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Understanding  the  behavior  of  soil  under  blast  loading  is  very  important  to 
engineers  in  mining,  tunneling,  and  military  construction.  Due  to  the  very  complex 
structure  of  a  soil  mass  it  is  very  difficult  to  describe  its  constitutive  relation,  especially 
when  it  has  different  water  contents  and  it  is  under  blast  loading  conditions.  New 
protective  system  designs  subjected  to  blast  loading  need  to  be  proved  its  validation  prior 
to  predict  effect  of  explosive  before  implementation.  Full-scale,  buried  explosive  tests  are 
costly.  Finite  element  simulations  play  a  significant  role  in  the  design  of  protective 
systems,  for  example  a  bottom  platform  of  lightweight  vehicles,  against  underground 
explosion. 

The  Perzyna  viscoplastic  cap  model  has  been  shown  to  be  a  valid  model  for  use  in 
the  simulations  of  dry  soil  behavior  under  both  static  and  dynamic  loading.  This  model  is 
a  dramatic  improvement  over  the  inviscid  cap  model  for  soil  behavior  under  high  strain 
rate  loading,  such  as  from  an  explosion.  However,  soil  should  be  modeled  as  a  three- 


phase  porous  media  to  accommodate  various  degrees  of  water  saturation.  This  is 
especially  true  for  the  soil  mass  surrounding  the  source  of  energy  release,  as  each  of  the 
three  phases  responds  differently  to  shock  loading.  To  improve  the  model  accuracy,  a 
revised  model  comprising  a  Gruneisen  equation  of  state  (EOS)  for  each  of  the  three 
phases  has  been  developed.  These  equations  of  state  for  solid,  water  and  air  have  been 
integrated  with  a  viscoplastic  cap  model  to  simulate  behaviors  of  soil  with  different 
degrees  of  water  saturation. 

These  EOS  models  as  well  as  the  viscoplastic  cap  model  are  implemented  into 
LS-DYNA  as  user-supplied  subroutines  for  numerical  simulation  of  six  explosive  tests  in 
dry  soil  as  well  as  in  saturated  soil.  The  shock  front  time  of  arrival,  the  air  pressure 
directly  above  the  buried  explosive,  and  the  ejecta  heights  predicted  by  the  revised  cap 
model  agree  fairly  well  with  the  experimental  data.  Four  elements  from  finite  element 
mash  are  selected  to  observe  three  phases  volume  fractions  change.  There  is  noticeable 
improvement  in  the  prediction  of  saturated  soil  behavior  than  dry  soil  behavior  under 
blast  loading.  It  is  concluded  that  the  revised  model  is  adequate  for  blast  loading  behavior 
simulations  for  soil  with  different  degrees  of  water  saturation. 
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CHAPTER  ONE  INTRODUCTION 


1 


1.1  BACKGROUND 

Many  commercial  and  military  endeavors,  such  as  defense,  construction, 
earthquake  prevention,  disaster  mitigation,  and  mining,  involve  soil  dynamics.  Soil 
behavior  under  blast  loading  have  been  studied  by  engineers  and  researchers  (Wang  and 
Lu  2003;  Tong  and  Tuan  2007;  Grujicic  et  al.  2008).  Soil  is  an  assemblage  of  individual 
particles,  rather  than  a  continuum,  that  soil  may  have  various  degrees  of  water  saturation. 
A  rapid  release  of  energy  from  a  buried  explosion  causes  a  sudden  rise  of  pressure  or  a 
shock  front  propagating  through  a  soil  medium,  it  is  very  challenging  to  accurately 
predict  soil  behavior  under  blast  loading.  Therefore,  to  date  common  practice  in  modeling 
soil  behavior  under  blast  loading  is  mainly  based  on  empirical  formulae  from  field  tests 
(Wang  et  al.  2004).  Since  conditions  varied  in  those  test  sites,  predictions  using  those 
empirical  formulae  scatter  significantly.  Discrepancy  at  the  same  scaled  distance  could 
be  more  than  two  orders  of  magnitude  between  dry  and  saturated  soils  (Drake  and  Little 
1983). 


Soil  is  composed  of  solid  particles  with  different  sizes  and  shapes  that  form  a 
skeleton  and  the  voids  are  filled  with  water  and  air.  The  soil  is  saturated  if  all  the  voids 
are  filled  with  water.  Otherwise,  the  soil  is  partially  saturated.  If  all  the  voids  are  filled 
with  air,  the  soil  is  said  to  be  dry.  It  is  a  common  practice  in  soil  mechanics  to  assume 
that  the  solid  particles  do  not  deform  and  the  water  phase  is  incompressible.  Hence, 
external  loading  is  supported  by  the  skeleton  and  the  water.  The  “effective  stress”  is  the 
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average  stress  on  a  plane  through  the  soil  mass,  rather  than  the  contact  stress  between  the 
solid  particles.  The  stress  on  the  water  and  the  air  is  called  “pore  pressure.”  The  principle 
of  effective  stress  was  first  recognized  by  Terzaghi  in  the  mid- 1920s  during  his  research 
in  soil  consolidation  (Budhu  2007).  Soils  cannot  sustain  significant  tension,  and  thus  the 
effective  stress  cannot  be  tensile.  Pore  water  pressures,  however,  may  be  positive  or 
negative  (i.e.,  suction).  For  unsaturated  soils,  the  effective  stress  (Bishop  1959)  is 
expressed  as 

a'=cr-Ua+z(ua-uw) . (1.1) 

where  a  is  the  total  stress,  ua  is  the  pore  air  pressure,  uw  is  the  pore  water  pressure,  and  x 
is  a  factor  depending  on  the  degree  of  saturation.  For  dry  soil,  x=0;  for  saturated  soil,  x=l 
(Loret  and  Khalili  2000;  Budhu  2007).  For  instance,  values  of  x  for  silts  are  shown  in  FIG. 
1-1. 


Degree  of  saturation  (%) 


FIG.  1-1  Values  of  for  a  silt  at  different  degrees  of  saturation 


3 

A  number  of  investigators  have  clearly  demonstrated  the  effective  stress 
hypothesis  under  static  and  quasi-static  loading  because  the  deformation  of  the  soil 
skeleton  depends  on  the  effective  stress  caused  by  the  structural  configuration  of  the  solid 
particles,  while  the  moisture  and  air  are  assumed  to  flow  through  the  skeleton  driven  by 
the  pore  pressure.  The  effective  stress  approach  becomes  invalid  under  shock  loading. 
This  is  due  to  the  fact  that  solid  particles  will  deform  under  shock  loading,  while  moisture 
and  air  are  trapped  in  soil  pores,  providing  additional  load  support. 

For  simulation  accuracy  in  finite  element  analysis  reasonable  constitutive  models 
for  the  involved  materials  are  critical.  Three  materials,  explosive,  air  and  soil,  are 
essential  to  define  an  underground  explosion.  The  constitutive  models  for  explosive  and 
air  have  been  reasonably  described  and  are  available  for  explosion  simulation  (Dobratz 
and  Crawford  1985,  “LS-DYNA”  2001),  but  soil  models  not  be  adequately  have 
implemented  into  finite  element  programs  for  explosion  simulations. 

1.2  OBJECTIVES  OF  THE  RESEARCH 

The  objective  of  this  research  is  to  develop  a  soil  model  developed  for  finite 
element  simulations  of  explosions  in  soil  with  various  degrees  of  saturation.  Equation  of 
state  (EOS)  models  are  developed  for  the  three  phases  of  the  soil  based  on  Kandaur’s 
concept  (Henrych,  1979).  These  EOS  models  are  integrated  with  the  viscoplastic  cap 
model  previously  developed  by  Tong  and  Tuan  (2007),  and  then  incorporated  into  LS- 
DYNA  as  user-defined  subroutines  for  soil  constitutive  relationship.  This  revised  cap 
model  is  then  validated  by  comparing  simulation  results  against  experimental  data. 
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Explosive  tests  conducted  by  Materials  Sciences  Corporation  (2006),  in  saturated  soil  as 


well  as  in  dry  soil,  were  used  to  validate  the  revised  cap  model. 

1.3  CONTENTS 

This  thesis  is  organized  as  follows. 

Chapter  One:  The  background,  objective  and  contents  of  this  study  are  described. 

Chapter  Two:  The  Perzyna  type  viscoplastic  cap  model  is  prepared  incorporating 
the  viscoplastic  cap  models  into  LS-DYNA  finite  element  code  as  user-defined  material 
models. 


Chapter  Three:  Two  formulations  of  equation  of  state  based  on  Kandaur 
conceptual  method  are  described.  An  equation  of  state  for  soil  is  established  and 
incorporated  the  equation  of  state  into  LS-DYNA  finite  element  code  as  user-defined 
equation  of  state  model. 

Chapter  Four:  The  models’  performance  is  evaluated  using  soil  viscoplastic  cap 
model  with  equation  of  state  in  finite  element  simulation  of  a  series  of  mine  explosion 
tests.  Four  elements  from  finite  element  mash  are  selected  to  observe  three  phases 
volume  fractions  change. 

Chapter  Five:  Conclusions  of  the  research  project  are  presented  as  well  as 
suggestions  and  recommendations  for  future  study. 
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CHAPTER  TWO  VISCOPLASTIC  CAP  MODELS 

2.1  INTRODUCTION 

Soil  has  generally  a  complex  structure  consisting  of  mineral  particles  which  form 
a  soil  skeleton.  The  interstertices  between  the  solid  particles  are  filled  with  air  and/or 
moisture.  In  general,  components  of  soil  are  solid,  water  and  air  and  called  three-phase 
soil.  A  soil  element  is  illustrated  in  FIG.  2-1. 


FIG.  2-1  A  schematic  element 

The  soil  skeleton  can  transmit  normal  stresses  and  shear  stresses  through  the  inter 
particle  contacts.  This  skeleton  of  grains  behaves  in  a  very  complex  manner  that  depends 


on  a  large  number  of  factors,  among  which  void  ratio,  partial  shape,  distribution  of  partial 
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size  and  confining  pressure  are  the  most  important  (Lade  2005).  When  the  pores  between 
the  solid  particles  are  filled  with  air,  the  soil  is  referred  to  dry  soil.  When  the  pores  are 
filled  with  water  containing  a  small  fraction  of  air  the  soil  is  called  saturated.  The  relative 
volume  fractions  of  the  three  constituent  materials  in  the  soil  are  generally  quantified  by 
the  porosity, a,  and  the  degree  of  saturation,  [j.  which  are  respectively  defined  as: 


a 


Li 

1/ 


(2.1) 


and 


Where  V  is  the  volume  of  void  (pores),  Vw  is  the  volume  of  water  and  V  is  the  total 
volume  of  the  three  phase  material. 


For  many  low  load  rate  processes,  the  overall  macroscopic  behavior  of  the  soil 
skeleton  may  be  defined  within  the  principles  of  continuum  mechanics,  making  it 
possible  to  simplify  the  modeling  and  apply  the  theories  and  methods  of  continuum 
mechanics. 


For  rapid  loading  conditions,  soil  models  incorporate  constitutive  models  of  the 


three  phases  all  required  to  define  soil  behavior. 


2.2  DEVELOPMENT  OF  SOIL  MODELS 
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2.2.1  SOIL  BEHAVIOR 

In  this  section,  different  characteristics  of  soil  behavior  are  discussed. 

(1)  Shear  strength  and  deformation  characteristics:  The  energy  applied  to  a  soil 
through  external  loads  may  both  overcome  the  frictional  resistance  between  the 
soil  particles  and  also  to  expand  the  soil  against  the  confining  pressure.  The  soil 
grains  are  highly  irregular  in  shape  and  have  to  be  lifted  over  one  another  for 
sliding  to  occur  (Das,  1983).  The  relationship  between  the  shear  strength  of  a  soil 
and  its  deformation  characteristics  depends  mainly  on  how  the  volume  changes 
during  the  shearing  process.  This  behavior  is  called  dilatancy.  An  increase  in 
volume,  or  expansion,  is  known  as  positive  dilation;  while  a  decrease  in  volume, 
or  contraction,  is  known  as  negative  dilation.  A  typical  curve  of  the  soil  dilatation 
under  shear  loading  is  shown  in  FIG.  2-2.  In  the  case  of  sands,  the  degree  of 
interlocking  between  particles  is  greater  when  the  soil  is  densely  packed.  An 
initial  expansion  or  dilation  is  necessary  in  order  for  the  soil  particles  to  more  past 
each  other.  Thus  the  shear  stress  will  first  rise  sharply  to  a  peak  value  at  a 
relatively  low  value  of  displacement,  with  a  corresponding  increase  in  volume.  At 
this  new  volume  the  interlocking  is  reduced  and  consequently,  as  the 
displacement  is  continued,  the  shear  stress  falls  back  and  finally  levels  off  at  an 
ultimate  residual  value  (Whitlow,  1995). 

(2)  Plasticity:  An  increase  in  applied  stress  usually  brings  about  some  irrecoverable 
deformation,  without  any  signs  of  cracking  or  disruption  (unloading,  see  FIG.  2-2). 
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Most  soils  only  have  a  very  small  elastic  region  and  show  plasticity  from  the 
onset  of  loading. 

(3)  Strain-hardening/softening:  After  an  initial  extension,  the  soil  behaves  as  if  it 
had  acquired  better  elastic  properties  and  a  higher  elastic  limit,  while  at  the  same 
time  it  had  lost  a  great  part  of  the  plastic  strain.  And  yield  surface  changes  with 
plastic  strain  development  during  loading  (Maugin,  1992).  Associated  with  strain¬ 
softening  behavior  is  the  tendency  of  dense  granular  and  overconsolidated  clays 
to  dilate  when  sheared  strain-hardening  is  associated  with  compaction  of  such 
materials  as  loose  sand  or  normally  consolidated  clays  experience  strain¬ 
hardening  (FIG.  2-2). 


t 


FIG.  2-  2  Response  of  soil  with  respect  to  shearing 
(Whitlow  1995) 
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(4)  High  Strain  Rate  Behaviors:  Soil  with  different  varying  water  contents  show 
different  behaviors  under  high  strain  rate  loading.  Test  data  using  a  Split 
Hopkinson  Pressure  Bar  (SHPB)  (Bragov  et  al.  2005;  Proud  et  al.  2007)  showed 
that  the  density  of  soil  and  the  shock  velocity  are  increased  with  moisture  content 
increasing.  A  schematic  SHPB  test  setup  is  shown  in  FIG.  2-3. 


FIG.  2-  3  Split  Hopkinson  Pressure  Bar  test  setup 


A  compressive  stress  pulse  is  induced  in  the  incident  bar  by  a  striker,  and  the 
incident  strain  ST ,  reflected  strain  £,,  and  transmitted  strain  £,r  in  the  bars  are 

IK’’  1 

measured.  The  stress-strain  relation  of  the  soil  specimen  and  the  strain  rate  can  be 
determined  from  the  elastic  modulus  of  the  bars  and  the  recorded  strain  data.  The 
confined  axial  stress-strain  curves  of  the  soil  specimens  from  SHPB  tests  at  three 
different  strain  rates  are  presented  in  FIG.  2-4. 
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FIG.  2-  4  Split  Hopkinson  Pressure  Bar  test  data 

(5)  Tensile  strength:  In  granular  media,  tensile  strength  is  a  result  of  various 
interparticle  physicochemical  forces  such  as  a.)  van  der  Waals  attraction,  b.) 
electrical  double  layer  repulsion  or  attraction,  c.)  cementation  due  to  solute 
precipitation,  d.)  capillary  stress  due  to  the  negative  pore  water  pressure,  and  e.) 
capillary  stress  due  to  the  surface  tension  of  liquid  (Lu  and  Likos,  2006).  The 
macroscopic  manifestation  of  these  forces  is  the  cohesive  behavior  shown  widely 
in  granular  media.  This  strength  can  play  an  important  role  in  stress  and  strain 
behavior.  Experimental  tensile  strength  results  for  the  silty  sand,  fine  sand,  and 
medium  sand  are  plotted  in  FIG  2-5,  FIG  2-6  and  FIG  2-7,  respectively,  as  a 
function  of  saturation  (Fu,  Wu  and  Tan,  2007). 
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Saturation  (%) 


FIG.  2-  5  Tensile  strength  curve  —  silt  sand 


Saturation  (%) 


FIG.  2-  6  Tensile  strength  curve  —  fine  sand 
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FIG.  2-  7  Tensile  strength  curve  —  medium  sand 


(6)  Effects  of  drainage  and  volume  change:  In  a  saturated  soil  mass,  an  increase  in 
applied  compressive  stress  of  creep  loading  causes  the  pore  water  pressure  to 
increase.  If  drainage  is  possible  an  outflow  of  water  then  takes  place  into 
surrounding  regions  where  the  pore  water  pressure  is  lower.  The  rate  of  outflow 
depends  on  the  permeability  of  the  soil,  in  gravels  and  sands  it  is  relatively  rapid, 
but  in  silts  and  clays  it  is  slow.  As  the  excess  pore  water  pressure  is  dissipated,  the 
applied  stress  is  transferred  from  pore  pressure  to  effective  stress.  Undrained 
conditions  occur  when  either  drainage  is  prevented  or  when  the  rate  of  application 
of  load  is  too  rapid  to  allow  significant  outflow  of  water.  The  deformation  of  an 
undrained  soil  mass  is  related  to  the  stiffness  of  both  the  pore  water  and  the  solids. 
When  loading  is  applied  slowly,  such  that  the  water  drains  away  without  any 
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increase  in  pore  water  pressure,  the  volume  will  decrease  and  stress-strain 

behavior  must  be  defined  in  terms  of  effective  stresses  (Whitlow,  1995). 

It  should  be  mentioned  that  there  are  also  other  characteristics  of  soil  behavior 
such  as  creep  and  temperature-dependency.  Those  aspects  are  not  discussed  here  because 
they  are  beyond  the  scopes  in  this  study. 

2.2.2  SOIL  MODELS 

The  mechanical  behavior  of  soils  may  be  modeled  at  many  levels.  Hooke's  law  of 
linear,  isotropic  elasticity  may  be  thought  of  as  the  simplest  available  stress-strain 
relationship,  but  this  is  generally  too  crude  to  capture  the  essential  features  of  soil 
behavior  (Brinkgreve,  2005).  On  the  other  hand,  a  large  number  of  constitutive  models 
have  been  proposed  by  several  researchers  to  describe  various  aspects  of  soil  behavior  in 
detail  (Lade  2005,  Prevost  and  Popescu  1996,  Chen  and  Baladi  1985).  Models  that  are 
appropriate  to  be  implemented  into  finite  element  programs  and  to  predict  the  soil 
behavior  for  geotechnical  engineering  applications  are  rather  limited. 

This  study  is  focused  on  a  limited  number  of  frequently  used  soil  models  that  can 
predict  the  soil  behavior  previously  described.  These  models  include  elastic  perfectly- 
plastic  soil  models,  hardening-plastic  soil  models,  elastic-viscoplastic  soil  models,  three- 
phase  soil  models,  viscoplasitc  soil  models,  SFG  (presented  by  Sheng,  Fredlund  and 
Gens)  unsaturated  soil  model  and  bounding  surface  plasticity  unsaturated  soil  models. 


2.2.2.1  ELASTIC  PERFECTLY-PLASTIC  SOIL  MODELS 


14 


The  classical  Mohr-Coulomb  model  is  often  used  to  describe  soil  behavior.  In  one 
dimension,  the  yield  “surface”  of  Mohr-Coulomb  mode  is  defined  by  a  linear  line 
between  shear  stress  r  and  normal  stress  cr  which  is  written  as 


/  =  r|  -  (c  -  a  tan  <f)  =  0 


.(2.4) 


where  the  constant  c  and  (f>  are  cohesion  and  internal  friction  angle.  But  in  three 
dimensions,  the  yield  surface  defined  by  Eq.  2.5,  is  much  more  complicated. 


•(2.5) 


/  =  ^  /j  sin  (j)  +  777  sin(<9  +  ^)  +  cos (0  +  ^)  sin  (j>  -  c  cos (j)  -  0 


where  1\  =  tr  <j(cr=  stress),  the  first  invariant  of  stress  tensor; 


J2  =  1/2  s  :  s  (s  =  stress  deviator),  the  second  invariant  of  deviatoric  stress  tensor; 


0  =  the  angle  of  similarity  and  defined  by  Eq.  2.6. 


cos  30  = 


3V3  J3 
2  Jln 


.(2.6) 


where  J3  =  det  Isl,  the  third  invariant  of  deviatoric  stress  tensor.  The  failure  surface  of 
Mohr-Coulomb  model  in  principal  stress  space,  which  is  hexagonal,  is  shown  in  FIG.  2-8. 


The  Mohr-Coulomb  model  is  certainly  still  the  best  known  model  for  an  isotropic 
pressure-sensitive  soil,  since  the  stress  at  failure  through  experimental  studies  agrees  well 
with  this  model  (Goldscheider,  1984).  The  model,  however,  is  not  mathematically 
convenient  due  to  the  presence  of  comers  or  singularities  (see  FIG.  2-8).  A  reasonable 
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smooth  generalization  of  the  Mohr-Coulomb  model  in  three  dimensional  situations  may 
be  defined  by  the  Drucker-Prager  model  (1952),  which  is  a  simple  cone  in  principal  stress 
space  as  shown  in  FIG.  2-9. 


FIG.  2-  8  Mohr-Coulomb  model  FIG.  2-  9  Drucker-Prager  model 

(Brinkgreve  2005)  (Brinkgreve  2005) 

Both  Mohr-Coulomb  model  and  Drucker-Prager  model  capture  soil  plasticity 
behavior  very  well  and  ensure  a  unique  solution.  However,  these  perfectly-plastic  soil 
models  have  inherent  limitations  and  shortcomings:  (1)  the  amount  of  dilatancy  predicted 
by  these  models  is  much  greater  than  observed  experimentally;  (2)  tests  indicate  a 
considerable  hysteresis  in  a  hydrostatic  loading-unloading  which  cannot  be  predicted 
using  the  same  elastic  bulk  modulus  of  loading  and  unloading  and  a  yield  surface  which 
does  not  cross  the  hydrostatic  loading  axis  (DiMaggio  and  Sandler,  1971);  (3)  strain 
softening  behavior  cannot  be  reproduced;  and  (4)  strain  rate  effect  is  not  considered. 


2.2.22  HARDENING-PLASTIC  SOIL  MODELS 
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Considering  strain  hardening/softening  behavior  of  soil,  it  is  logical  to  utilize  the 
classical  theories  of  work-hardening  plasticity  developed  for  metals.  Drucker  et  al.  (1957) 
first  suggested  that  soil  might  be  modeled  as  an  elastic -plastic  work-hardening  material. 
They  proposed  that  successive  yield  surfaces  might  resemble  extended  Drucker-Prager 
cone  with  convex  end  spherical  caps  as  shown  in  FIG.  2-10  (Chen  and  Baladi  1985).  As 
the  soil  strain  hardens,  both  the  cone  and  the  end  cap  expand.  This  concept  of  cap 
envelope  was  a  major  step  forward  toward  a  more  realistic  representation  of  soil  behavior. 
Based  on  this  concept,  numerous  work-hardening  soil  models  have  been  developed. 
Generally  these  models  can  be  classified  into  two  groups:  modified  Cam-clay  model  and 
generalized  cap  model. 


FIG.  2-  3 


A  Drucker-Prager  type  of  strain-hardening  cap  model 
(Chen  and  Baladi  1985) 
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The  modified  Cam-clay  model  was  developed  at  Cambridge  University  by 
Roscoe  et  al.  (1963).  This  model  is  based  on  critical  state  theory  and  originally  meant  to 
simulate  the  behavior  of  near-normally  consolidated  clays  under  triaxial  compression  test 
conditions.  The  fundamental  concept  of  this  model  is  that  there  exists  a  unique  failure 
surface  that  defines  failure  of  a  soil  irrespective  of  the  history  of  loading  or  stress  paths 
followed.  The  yield  surface  is  assumed  to  be  an  ellipse  and  may  be  expanded  with  the 
increase  of  volumetric  strain,  as  shown  in  FIG.  2-11  in  the  stress  space  of  I\  ~  C/2- 


FIG.  2-  4  Modified  Cam-Clay  model 

By  taking  Drucker-Prager  yield  surface  as  the  critical  state,  the  Cam-clay  models 
can  not  only  predict  the  failure  behaviors  very  well,  but  also  describe  the  nonlinear  and 
stress-path  dependent  behaviors  prior  to  failure  accurately,  especially  for  clay-type  soils. 
This  model,  however,  still  has  some  disadvantages  (DiMaggio  and  Sandler,  1971):  (1)  the 
discontinuous  slope  at  the  intersection  with  the  I\  axis  predicts  behavior  that  does  not 
seem  to  be  supported  by  experiments;  (2)  points  on  the  yield  surface  above  the  critical 
state  line  do  not  satisfy  Drucker’s  postulate  of  stability. 
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The  generalized  cap  model  was  first  proposed  by  DiMaggio  and  Sandler  (1971) 
based  on  the  concept  of  Drucker  et  al  (1957).  The  yield  function  consists  of  a  perfectly- 
plastic  (failure)  portion  fitted  to  a  strain-hardening  elliptical  cap  as  shown  in  FIG.  2-12. 
The  movement  of  the  cap  surface  is  controlled  by  the  increase  or  decrease  of  the  plastic 
volumetric  strain,  strain-hardening  can  therefore  be  reversed.  It  is  this  mechanism  that 
leads  to  an  effective  control  of  dilatancy,  which  can  be  kept  quite  small  as  required  for 
many  soils.  Moreover,  the  functional  forms  for  both  the  perfectly-plastic  and  strain¬ 
hardening  portions  may  be  quite  general  and  allow  for  the  fitting  of  a  wide  range  of 
material  properties.  With  the  associated  flow  rule,  this  model  may  satisfy  all  the 
theoretical  requirements  of  stability,  uniqueness  and  continuity.  The  agreement  between 
model  behaviors  and  static  experimental  results  for  granular  soils  are  reasonably  good. 
As  for  the  disadvantages  for  this  model,  one  is  the  corners  at  the  intersection  of  the  yield 
curves  which  may  cause  mathematical  problems.  If  the  stress  status  happens  to  fall  in  the 
comers,  the  consistency  condition  may  not  be  fulfilled. 


FIG.  2-  5  Yield  surface  of  generalized  cap  model 
(DiMaggio  and  Sandler  1979) 
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Numerous  formulations  have  been  proposed  in  the  literature  to  improve  the 


capacity  of  this  model.  Sandler  et  al  (1976)  introduced  a  more  generalized  cap  model 
with  different  modulus  on  loading  and  unloading.  Later  Sandler  and  Rubin  (1979) 
devised  a  cap  model  algorithm  which  permitted  to  obtain  flexibility  with  respect  to 
changes  in  functional  forms  and  parameters.  Simo  et  al  (1988)  proposed  a  new  algorithm 
in  which  special  attention  was  paid  to  the  singular  corner  regions  at  the  intersection  of  the 
yield  surfaces  in  order  to  be  consistent  with  the  notion  of  the  close-point  projection 
method.  Various  smooth  cap  models  were  also  proposed  to  eliminate  the  numerical 
problem  at  corners  (Rubin  1991,  Schwer  and  Murry  1994).  The  continuous  surface  cap 
model  developed  by  Rubin  (1991)  is  shown  in  FIG.  2-13. 


FIG.  2-  6  Stress  Space  View  of  Continuous  Surface  cap  model 


(Rubin  1991) 
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2.2.23  THREE-PHASE  SOIL  MODELS 

In  the  early  1980s,  the  development  of  constitutive  equations  for  saturated  soils 
required  three  main  components:  the  concept  of  effective  stress;  elastic -plastic  models 
based  on  effective  stress  able  to  describe  the  behavior  of  drained  soil  under  complex 
loading  paths  and  finally,  the  theory  of  mixtures  for  a  solid  skeleton  and  fluid.  Loret  and 
Khalili  (2000)  developed  a  three-phase  model  for  unsaturated  soils  which  despite  the 
similarity  of  the  framework  presented.  There  are  numerous  differences  between  the 
saturated  and  unsaturated  soils.  For  saturated  soils,  the  concept  effective  stress  developed 
by  Terzaghi  is  seldom  questioned  in  its  role  describing  the  plastic  behavior  of  saturated 
soils.  The  situation  is  far  from  being  identical  for  unsaturated  soils,  which  are  three-phase 
materials.  Bishop  and  Blight  (1963)  provided  experimental  evidence  supporting  the 
validity  of  Bishop’s  stress,  they  observed  that  the  volumetric  characteristics  and  shear- 
strength  do  not  change  if  the  effective  stress  for  the  individual  components  vary  in  such  a 
way  as  to  keep  the  net  stress  and  suction  constant.  However,  Jennings  and  Burland  (1962) 
questioned  the  validity  of  the  principle  of  effective  stress  in  the  deformation  behavior  of 
unsaturated  soils  arguing  that  it  cannot  explain  the  collapse  phenomenon  observed  during 
wetting.  On  the  other  hand,  they  agree  that  the  shear-strength  depends  on  an  effective 
stress  of  the  form  proposed  by  Bishop.  Burland  (1965)  further  resorted  to  arguments  of 
theoretical  nature  reasoning  against  adding  a  macroscopic  quantity,  the  net  stress,  and  a 
microscopic  quantity,  the  suction.  Although  this  was  not  checked,  their  arguments  have 
been  widely  accepted  and  several  researchers  have  concluded  that  the  behavior 


description  of  unsaturated  soils  should  consider  net  stress  and  suction  as  two  independent 
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variables.  Effective  stress  is  the  key  quantity  that  governs  the  material  behavior,  both  in 
the  elastic  and  plastic  regimes. 


The  third  invariant  of  the  effective  stress  is  not  accounted  for  and  the  cross-sections  along 
deviatoric  planes  are  circular,  FIG2-14, 


/  =  f(p,q,Pc)  =  -£r=  +  P~  Pc  . (2-9) 

M-p 

where  M  is  the  slope  of  the  critical  state  line;  it  is  assumed  to  be  a  material  parameter, 
independent  in  particular  of  suction.  The  size  of  the  yield  surface  can  be  measured  by  the 
pre-consolidation  stress  Pc . 


FIG.  2-  7  Yield  Surface  of  the  Modified  Cam-Clay  Model  in  terms  of  the  Effective 


Mean-Stress  and  Shear  Stress  (Foret  and  Khalili  2000) 
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However,  the  three-phase  soil  model  has  the  limitations  and  shortcomings:  (1) 


experimental  results  are  not  ready  available  to  justify  the  use  of  the  complex  model;  and 
(2)  the  identification  of  material  coefficients  require  the  use  of  experiments  ,for  example, 
the  soil-water  characteristic  curve  is  incorporated  into  the  model. 

2.2.2.4  VISCOPLASTIC  SOIL  MODELS 

Viscoplastic  models  are  expansion  of  plastic  models  which  include  rate  effects.  If 
the  yield  surfaces  in  viscoplastic  model  are  taken  the  same  as  those  in  plastic  model,  a 
simple  explanation  of  the  difference  between  viscoplastic  and  plastic  solution  may  be 

shown  in  FIG.  2-15.  In  the  stress  space  of  7j  ~  V. J2 ,  the  plastic  solution  ( CFn+l )  must  be  lie 

on  one  of  the  yield  surfaces,  this  is  violated.  The  viscoplastic  solution  ( <Jn+l )  may  be 
outside  of  the  yield  surface  due  to  the  rate  effect.  From  the  point  view  of  numerical 
calculation,  plasticity  may  be  considered  as  a  special  case  of  viscoplasticity  when  strain 
rate  is  low  enough  to  be  neglected.  Viscoplasicity  is  an  improvement  of  plasticity  in  its 
ability  to  predict  the  soil  behavior  over  a  wide  range  of  loading  rate. 


FIG.  2-  8  Viscoplasticity  vs.  plasticity 
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The  most  well-known  formulation  of  viscoplasticity  is  based  on  Perzyna’s  theory 


(1966),  wherein  viscous  behavior  is  modeled  with  a  time-rate  flow  rule.  The  flow  rule  is 
assumed  to  be  associative  such  that  the  viscoplastic  potential  is  identical  or  at  least 
proportional  to  the  yield  surface  (Katona  1984,  Chen  and  Baladi  1985,  Simo  et  al.  1988). 
Perzyna  (1966)  pointed  out  that  the  models  with  rate-dependent  elastic  response  (i.e. 
viscoelastic  models  such  as  Murayama  and  Shibata  (1964))  are  mathematically  very 
complicated.  In  addition,  rate-independent  elastic  response  models,  because  of  their 
relative  mathematical  simplicity  and  their  similarities  with  the  inviscid  theory  of 
plasticity,  may  be  more  appropriate  for  practical  engineering  applications  (Perzyna  1966, 
Swift  1975).  Also,  viscous  effects  appear  to  be  more  evident  in  the  plastic  range  for  most 
clay.  Models  which  explicitly  incorporate  time  into  the  stress-strain  relations  (Matsui  and 
Abe  1985,  Sekiguchi  1984)  have  the  drawback  of  predicting  time-dependent 
deformations  under  zero  effective  stress  condition.  Also,  it  is  difficult  to  determine  the 
correct  value  of  the  material  time  parameter  if  the  stress  history  is  not  known.  Dafalias 
(1982),  from  microscopic  and  macroscopic  observations  of  the  structure  and  response  of 
clays,  has  concluded  the  constitutive  relations  can  best  be  obtained  by  considering  the 
plastic  strains  as  a  combination  of  rate-dependent  and  non-rate  dependent  components 
(elastoplastic  viscoplastic  models  such  as  those  of  Kaliakin  (1985)  and  Broja  and 
Kavazanjian  (1985)).  However,  there  is  a  difficulty  in  this  formulation.  Beyond  a  cetain 
strain  rate,  further  increases  do  not  affect  the  stress-strain  relationship  (Dafalias,  1982). 
Effects  of  very  high  strain  rates  cannot  therefore  be  predicted  (Whitman  1957, 
Richardson  and  Whitman  1963,  Adachi,  Mimura  and  Oka  1985).  Although  the 
viscoplastic  model  of  the  Perzyna  type  has  been  validated  through  simple  dynamic  tests, 
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little  research  work  has  attempted  to  apply  this  model  to  simulate  the  soil  behavior  under 
an  extremely  high  strain-rate  loading  such  as  explosions. 

Another  popular  format  of  viscoplasticity  is  based  on  Duvant-Lions’  theory 
(1972),  wherein  the  viscoplastic  solution  is  simply  constructed  through  the  relavent 
plastic  solution.  The  biggest  advantage  of  the  Duvant-Lions’  model  is  its  ease  in 
numerical  implementation,  only  a  simple  stress  update  loop  is  needed  to  add  into  existing 
plasticity  algorithms.  Also  deterioration  from  viscoplastic  solution  to  plastic  solution 
under  a  low  strain  rate  is  exactly  guaranteed  (Simo  et  al  1987).  Compared  to  the  Perzyna 
type,  the  viscoplastic  model  of  the  Duvant-Lions’  type  has  not  been  well  validated 
experimently. 

2.2.2.5  SFG  UNSATURATED  SOIL  MODEL 

Since  the  pioneering  work  of  Alonso  et  al.  (1990),  a  number  of  elastoplastic 
constitutive  models  have  been  developed  for  modeling  the  behavior  of  unsaturated  soils 
(Gens  1996;  Jommi  2000;  and  Gens  etal.  2006).  Early  models  dealt  only  with  the  stress- 
suction- strain  relationships  of  unsaturated  soils  (Kohgo  et  al.  1993;  Wheeler  and 
Sivakumar  1995;  Bolzon  et  al.  1996;  Cui  and  Delage  1996;  Loret  and  Khalili  2002). 
These  models  are  based  on  the  same  basic  assumptions  and  largely  fall  in  the  same 
framework  as  Alonso  et  al.  (1990),  although  different  constitutive  equations  and  different 
stress  variables  are  used.  The  model  by  Alonso  et  al.  (1990),  generally  referred  to  as 
Barcelona  Basic  Model,  remains  as  one  of  the  fundamental  models  for  unsaturated  soils. 
More  recent  models  have  incorporated  suction-saturation  relationships  with  hysteresis 
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(Wheeler  1996;  Dangla  et  al.  1997;  Vaunat  et  al.  2000;  Gallipoli  et  al.  2003;  Wheeler  et 
al.  2003;  Sheng  et  al.  2004;  Santagiuliana  and  Schrefler  2006;  Sun  et  al.  2007). 

Usually,  unsaturated  soil  models  use  a  load-collapse  yield  surface  to  define  the 
variation  of  the  apparent  pre-consolidation  stress  along  the  soil  suction  axis.  The  apparent 
pre-consolidation  stress  is  usually  assumed  to  increase  with  increasing  suction.  Under 
such  a  framework,  these  models  are  able  to  reproduce  some  basic  features  of  unsaturated 
soil  behavior,  for  example,  the  volume  change  upon  wetting  and  the  increase  of  shear 
strength  with  suction.  However,  some  basic  questions,  like  how  yield  stress  changes  with 
soil  suction,  have  not  been  fully  answered.  The  SFG  model  presents  a  volumetric 
behavior  model  for  independent  changes  of  mean  net  stress  and  suction.  Based  on  this 
volumetric  relationship,  the  change  of  the  yield  stress  with  suction  and  the  hardening 
laws  that  govern  the  evolution  of  the  yield  surface  are  derived.  Recent  developments  have 
included  combining  both  stress-strain  and  suction-saturation  relations  of  unsaturated  soils 
are  also  incorporated  into  this  model. 

The  SFG  model  is  expressed  in  the  plane  mean  net  stress  p  and  suction  s  with 

P  =  P~ua . (2-10) 

s  =  ua-uw . (2.11) 

where  ua  is  the  pore  air  pressure  and  uw  is  pore  water  pressure. 


A  yield  surface  can  be  expressed  as, 
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.(2.12) 


where  Ssa  is  the  saturation  suction  and  pyi0  is  the  new  yield  stress  at  zero  suction.  If  pyt0 
is  known,  Eq.  (2-10)  can  be  used  to  find  the  new  yield  surface,  p .  Alternatively,  if  the 


new  yield  stress  at  a  given  suction  is  known,  Eq.  (2-10)  can  be  used  to  find  pyi0  . 


The  new  yield  surface  pyn  for  pyt0  =500  kPa  is  shown  in  FIG.  2-16.  The  yield 

stress  along  the  new  yield  surface  does  not  monotonically  decrease  with  increasing 
suction. 


FIG.  2-  9  Initial  Yield  Surface  for  a  Soil  That  was  Consolidated  to  300  kPa  at  Zero 
Suction  and  Its  Evolution  When  the  Soil  is  then  Loaded  at  Different  Suction  Levels 


(Sheng,  D.,  Fredlund,  D.G.  and  Gens,  A.  2008) 
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However,  the  SFG  soil  model  has  limitation  and  shortcomings:  (1)  there  are  not 
enough  experimental  data  at  present  to  define  precisely  the  suction-dependence  of 
material  parameters  and  (2)  the  yield  surface  function  is  different  with  different  suction 
level. 

2.2.2.6  BOUNDING  SURFACE  PUASTICITY  UNSATURATED  SOIU  MODEU 

Bounding  surface  plasticity  was  first  developed  for  metals  by  Dafalias  and  Popov 
(1976),  and  later  applied  to  clays  by  Dafalias  and  Herrmann  (1982),  to  pavement  base 
materials  by  McVay  and  Taesiri  (1985),  and  to  sands  by  Hashigushi  and  Ueno  (1977), 
Aboim  and  Roth  (1982)  and  Bardet  (1985).  By  broadening  the  scope  of  conventional 
plasticity  theory,  bounding  surface  plasticity  provides  a  flexible  theoretical  framework  to 
model  the  cyclic  behavior  of  engineering  materials. 

The  bounding  surface  plasticity  soil  model  represents  the  macroscopic  behavior  of 
soil  in  terms  of  strain  and  effective  stress  and  is  useful  by  the  solution  of  boundary  value 
problems  with  finite  element  or  finite  difference  methods.  The  advantages  of  bounding 
surface  plasticity  theory  over  conventional  plasticity  have  investigated  for  cyclic  as  well 
as  monotonic  loadings.  The  existence  and  direction  of  the  irreversible  strain  increment, 
which  requires  non-associative  the  flow  rule  to  be  realistically  simulated  at  the  failure 
state,  can  be  defined  by  only  one  surface  in  bounding  surface  plasticity. 

Russell  and  Khalili  (2005)  developed  an  unsaturated  soil  model  using  bounding 
surface  plasticity.  However,  this  model  took  particle  crushing  into  account,  making  it 
complex  and  difficult  to  apply  for  ordinary  cases  in  soil  mechanics.  Wong,  Morvan  and 
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Branque  (2009)  developed  a  new  bounding  surface  plasticity  model  for  unsaturated  soils 
with  a  small  number  of  parameters  based  on  Bardet’s  model  (Bardet,  1985).  In  this  model, 
the  bounding  surface  itself  evolves  through  a  hardening  mechanism  that  depends  on  the 
accumulated  plastic  strains. 

The  bounding  surface  of  this  model  is  elliptical  in  the  plane  effective  mean-stress 
p'and  shear- stress  q  with 


P'  =  \(a  l'+<72  +Cr3) 


—  ^1  ®  2 . 

The  bounding  surface  can  be  defined  as  (2-12),  FIG.  2-17, 


.(2.13) 

.(2.14) 
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where, 

P=A . (2-16) 

q'  =  pcM K  A  ^ . (2.17) 


x  = 


q 

Mp'  +  q  o 


(2.18) 


l  +  (p-l)Vl  +  ^V(p-2) 

1  +  (p  - 1)2  x2 


(2.19) 


Mk  is  the  slope  of  the  saturated  soil  critical  state  line  (CSL).  The  size  of  the  yield  surface 
can  be  measured  by  the  hardening  variable  AK.  Mn  and  An  are  assumed  to  be  a  material 


parameter,  independent  in  particular  of  suction  s.  p  is  a  material  parameter. 
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FIG.  2-  10  Yield  Surface  of  the  Modified  Unsaturated  Soil  Model  in  terms  of  the 
Effective  Mean-Stress  and  Shear  Stress 
(Wong,  Morvan  and  Branque  2009) 

The  bounding  surface  plasticity  soil  model  has  the  limitations  and  shortcomings: 
(1)  there  are  not  enough  experimental  data  at  present  to  precisely  define  the  suction- 
dependence  of  material  parameters  and  (2)  the  water  retention  curve  in  general  does  not 
define  an  objective  relation  between  degree  of  saturation  and  suction. 

2.3  VISCOPLASTIC  CAP  MODELS 

Viscoplasticity  is  defined  as  a  rate-dependent  (as  opposed  to  inviscid  means  rate 
independent)  plasticity  model  and  may  be  applied  to  the  soil  constitutive  laws  to  account 
for  the  strain  rate  effect.  A  variety  of  viscoplastic  formulations  for  soils  have  been 
proposed  in  the  literature.  The  formulation  of  viscoplasticity  based  on  Perzyna’s  theory 
(1966)  is  the  most  well-known  formulation,  where  viscous  behavior  is  modeled  with  a 
time-rate  flow  rule.  The  flow  rule  is  assumed  to  be  associative  such  that  the  viscoplastic 


30 

potential  is  identical  or  at  least  proportional  to  the  yield  surface  (Katona  1984,  Chen  and 
Baladi  1985,  Simo  et  al.  1988).  After  transition  into  rate-independent  plasticity,  this 
identity  becomes  essential  although  it  has  no  great  significance  in  viscoplasticity.  The 
viscoplastic  formulation  has  the  following  advantages:  (1)  the  generality  of  the  viscous 
flow  rule  offers  the  capability  of  simulating  time-dependent  material  behavior  over  a 
wide  range  of  loading;  and  (2)  the  extension  of  an  inviscid  cap  model  for  viscoplasticity 
is  relatively  straightforward  (Tong,  2005). 

Another  viscoplastic  formulation  of  the  Duvant-Lions  type  has  been  advocated  by 
Simo  et  al  (1988),  Simo  and  Govindjee  (1991)  and  Simo  and  Hughes  (1998).  The  viscous 
behavior  is  constructed  directly  based  on  the  difference  between  solutions  for  inviscid 
and  the  viscoplastic  foumulations.  The  main  advantage  is  its  ease  in  numerical 
implementation,  only  a  stress  update  needs  to  be  added  in  an  inviscid  formulation  in 
order  to  obtain  the  corresponding  viscoplastic  solution. 

The  viscoplastic  cap  model  is  an  effective  material  model  to  simulate  soil 
behavior  under  high  strain  rate  loading.  Tong  (2005)  applied  viscoplastic  cap  model  in 
LS-DYNA  to  simulate  a  series  of  explosions  in  soil.  Comparisons  with  experimental 
results,  the  simulations  of  soil  ejecta,  crater  and  explosive  clouds  from  landmine- 
explosion  tests  were  reasonably  good. 

In  this  chapter,  two  type  of  viscoplastic  cap  models  are  proposed  based  on 
Perzyna’s  theory  and  Duvant-Lions’  theory.  The  plastic  yield  functions  are  patterned  on 
the  generalized  two-invariant  cap  model.  Numerical  algorithm  is  presented.  The 
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performance  of  viscoplastic  cap  model  is  examined  using  a  hypothetical  uniaxial  strain 
test  and  compared  against  experimental  data  under  rapid  loading. 

In  the  viscoplasticity  model,  the  total  strain  rate  vector  &  is  decomposed  into  an 
elastic  component  M  and  a  viscoplastic  component  Mp 


&!’ . (2.20) 

The  elastic  component  is  expressed  as 

<&=  CaS . (2.21) 


whereat  =  the  stress  rate  vector  and  C  =  an  elastic  constitutive  matrix. 

For  the  viscoplastic  component,  it  is  different  with  the  different  type. 


2.3.1  THE  PERZYNA  TYPE  VISCOPLASTIC  CAP  MODEL 


The  viscoplastic  strain  rate  vector  is  assumed  to  be  delayed  with  time  and  is 
expressed  as  follows  when  assuming  associated  flow  rule: 


W  =r,<0(f)>^- . (2.22) 

da 

where  rj  =  a  material  constant  called  fluidity  parameter; 

x  +  lx 

The  notion  <  >  refers  to  the  ramp  function  defined  by  <  x  >=  — — — ; 

/  =  plastic  yield  function; 


0(f)  =  dimensionless  viscous  flow  function  and  commonly  expressed  in  the  form  of 
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where  N  =  an  exponent;  and/o  =  a  normalizing  constant  with  the  same  units  as/. 


2.3.1. 1  STATIC  YIELD  FUNCTIONS 


The  plastic  yield  function  /  is  patterned  in  the  inviscid  cap  model  (DiMaggio  and 
Sandler  1971,  Sandler  and  Rubin  1979,  Simo  et  al  1986)  which  is  formulated  in  terms  of 
the  first  stress  invariant  1\  and  the  second  deviator  stress  invariant  J2  as  shown  in  FIG.  2- 
18  (Tong,  2005).  The  static  yield  surface  is  divided  into  three  regions: 

(a)  when  /  >  L,  the  cap  surface  region  /  =  -Fc(I1,k)  =  0 

(b)  when  L  >  I\  >  -T ,  the  failure  surface  region  /  =  -  Fe  (/  )  =  0 

(c)  when  I\  <  -T ,  the  tension  cutoff  region  /  =  /  -  ( -T )  =  0 


FIG.  2-11  Static  yield  surface  for  cap  model  (Tong,  2005) 
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(a)  Cap  surface  portion:  the  cap  surface  is  a  hardening  surface  in  the  shape  of  an  ellipse 
quadrant  in  the  stress  space  of  7i  and  J2.  It  is  generally  defined  by 

f{lx,JJ~2,k)  =  -  FJ Ivk )  =  V77 ~^{x{k)- L(k)f  -(/1-  L{k)f  =0 . (2.24) 

where  h  =  tr  a  =  on  +  G22+  CJ33;  J2  =  1/2  s  :  s  =1/2  (sf l  2  +  $33 )  (s  =  stress  deviator); 

Fc(Ij,k )  =  the  loading  function  for  cap  envelope; 

R  =  a  material  parameter; 

k  is  a  hardening  parameter  related  to  the  actual  viscoplastic  volumetric  change 

s?(=ttsv  =  6'*  +  eZ+e$y. 


e?(Xm  =  w{l- exp[-  D(x(k)~  X  . )]; . (2.25) 

where  X(k )  defines  the  intersection  of  the  cap  with  the  7)  -  axis  and  hence  is  given  by 

X(k)  =  k  +  R-Fe(k) . (2.26) 

where  Fe(k )  =the  loading  function. 


L(k)  is  the  value  of  I\  at  the  location  of  the  start  of  cap  and  is  defined  by 

k)\k  if  k>  0 . (2.21) 

[0  if  k<0 

The  cap  surface  may  be  expressed  alternatively  (Katona  1984)  as 


/  (Ix,J2,k) 


(/|  ~  k)2  |  j  (L-X) 
R2  2  R 


(2.28) 


(b)  Failure  surface  portion:  the  failure  surface  is  a  non-hardening,  modified  Drucker- 


Prager  form  with  a  yield  function  defined  as 


/(vV^)=VA-^(A)  =  V^-(«-yex  p(-/?/,)+tf /,)  =  () 
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(2.29) 


in  which  a,  fJ.  /and  0  are  material  parameters. 

(c)  Tension  cutoff  portion:  the  tension  cutoff  surface  is  defined  by 

f  (I i)  -  1 1  ~  (~T)  . (2.30) 

where  -T  =  tension  cutoff  value. 

2.3.1.2  SOLUTION  ALGORITHMS 

The  strain  rate  Eq.  (2.20)  and  (2.21)  are  integrated  over  a  time  step,  At,  from  t  to 
t+At,  to  yield  the  incremental  strains  and  stresses: 


As  =  Ase  +  Asvp . (2.31) 

A<t  =  CAse  =  c(As-Asvp) . (2.32) 


where  As  =  the  total  incremental  strain  vector; 

A se  =  the  elastic  viscoplastic  incremental  strain  vector; 

A s'p  =  the  viscoplastic  incremental  strain  vector; 

Act  =  the  incremental  stress  vector. 

Based  on  the  Euler  method,  the  viscoplastic  incremental  strain  vector  As'7’  can  be 
approximated  as 

Ae  =  [(i-xW+x4L,Ut  . (2.33) 

in  which  %  is  an  adjustable  integration  parameter,  0  <  %  <  1.  The  integration  scheme  is 
explicit  if  %  =  0  and  fully  implicit  if  j  =  1.  This  solution  algorithm  is  conditionally  stable 
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when  %<  0.5  and  unconditionally  stable  when  %  >  0.5.  The  fully  implicit  integration 
scheme,  %=  1,  is  used  here  in  the  numerical  algorithm  just  for  simplification. 

In  the  full  implicit  integration  scheme,  the  viscoplastic  flow  (Eq.  2.33)  is  only 
determined  by  the  gradient  of  the  yield  surface  at  time  t+At.  Thus,  As'1’  may  be  rewritten 

as 

A e"p  =  A&pAt  =  rj  <  </>(f)>  At^- . (2.34) 

dcr 

If  a  plastic  multiplier  AA  is  introduced  such  that 

AT  =  ri  <  >  At . (2.35) 

then  Eq.  (2.34)  may  be  rewritten  as 

Ae'p  =AA^~ . (2.36) 

dc 

This  viscoplastic  problem  can  be  solved  under  the  condition  that  the  residual  p,  defined 
in  Eq.  (2.37),  is  reduced  to  zero  during  a  local  iteration. 


P  =  — - ^(/)->  0 . (2.37) 

rjAt 

Substituting  Eq.  (2.36)  into  Eq.  (2.32)  yields 

Act  =  C:(Ae-AA^-) . (2.38) 

dcr 


To  compute  AA,  a  local  Newton-Raphson  iteration  process  is  applied.  Note  that  the 
yield  function  takes  the  general  form  f  -  f(cr,k)  .  Differentiating  Eq.  (2.38)  during 


iteration  i  gives 
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r)f  d2 

da  =  C:  (Ss  -SX  —  -  AX<n  — — XrSa  -  AX 
8a  da1 


d'f  a  no  d2  f 


SX) 


.(2.39) 


dadX 

where  8cr,  8s  and  8X  are  the  iterative  improvements  of  Zlcr.  As  and  AX,  respectively, 
within  the  local  iteration  process. 


Eq.  (2.39)  may  be  expressed  alternatively  as 

Sa  =  U: 


Ss-(—  +  AX 


(0  8  f  )SX 


da  dad/ 1 

with  a  pseudo-elastic  stiffness  matrix  H 


.(2.40) 


H 


'■-l  ,  A  1  (0  ^  / 


-1-1 


C  +  A2 


da 


.(2.41) 


By  differentiation  of  Eq.  (2.38),  the  Newton-Raphson  process  at  iteration  i  is  expressed  as 

. (2.42) 


Pin  = 


_ d£ 

ijA  t  SX  y 


(  3 


8X- 


V  da  j 


8a 


Substitution  of  Eq.  (2.40)  into  Eq.  (2.42)  yields 


8X  =  - 
C 


with 


\da  j 


H  Ss  +  p 


w 


.(2.43) 


z= 


'df) 

T 

H 

\V+AA"dA  1 

,8a  j 

da  dadX 

+  ■ 


1  d<i> 


i]At  dX 


.(2.44) 


If  a  local  iteration  is  applied,  the  iterative  strain  increment  8s  will  turn  to  a  fixed 
total  strain  increment  As  during  a  global  iteration.  The  iterative  algorithm  for  the 
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viscoplastic  cap  model  of  the  Perzyna  type,  where  the  subscript  n  and  n+ 1  denotes  the 
solutions  at  time  t  and  t+At  respectively,  is  outlined  in  Table  2-1. 


Table  2-  1  Numerical  algorithm  for  the  Perzyna’ s  viscoplastic  model 


DATA  INPUT 

:  °„’kn,Ae 

Trial  stresses 

:  cj™  =an+CAe  ,kn 

If  f{<T%,km)<  0- 

elastic  a „+1  =  ,  U,+  i  =  kn  RETURN 

If  /(<?rf„)<0, 

viscoplastic 

(a)  define  the  initial  iteration  value 

AA(0)  =0  , 


cC'i  =  o’,,  +  C 


pm=t(°ZK)- 

(b)  do  local  iteration  i 

H 


,(0)  _ 

AT(0) 

r/At 


Ae-A^V- 

d<j 


(T1  +A2<,')  8  f 


da' 


T 

:  H : 

i  +  AA“> 

a2/ " 

va<7y 

dcr 

dadX 

4  = 


AX(i+l)  =  ax 


+  • 


1 


d(j) 
r/At  dX 


(0  .  P 


0) 


dX 


h.U+\)  _  1.(0  ,  ijL(i+l) 

%+l  —  Kn+1  +  mn+l 


_(/+!) 


P 


n+ 1  —  +  C  I 

=4 


As-AXM)dt- 

da 


O+D 


-—0+1)  7,  0+1) 

■  n+1  ’^n+l  , 


AT 


(i+D 


?7Ar 


go  back  and  continue  until  I//' 1  /'l  <  5 
RETURN 


OUTPUT 


^/i+l  ’  kn+ 1  >  ^n+1 


For  the  tension  cutoff  region,  a  different  algorithm  is  needed  since  the  yield 
surface  for  tension  cutoff  is  independent  of  J1.  This  condition  is  uncommon  for  ordinary 
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soil  tests  but  must  be  accommodated  for  explosions  in  soil.  Since  very  little  test  data  is 
available,  the  following  assumptions  are  made:  (1)  the  viscous  flow  parameter  under 
tension,  qv  may  be  the  same  as  or  different  from  that  under  compression;  (2)  the 

viscoplastic  solution  <T;+/V  is  presumed  to  be  between  the  elastic  trial  stress  and  the 

inviscid  solution  at+At ,  and  the  simplest  scheme  is  to  assume  that  <Jt  h/V  is  on  the  straight 

line  from  to  <r,+A; ,  as  shown  in  FIG.  2-19.  The  treatment  for  tension  cutoff  is  thus 
proposed  as  follows: 

(1)  if  Ci,  <-T  and  <  Fe (-T) ,  then 

T  —  -r/r At  j trial  ,(]_  ~>lr* <V_T\.  [J  _  /  J  trial 

1  l,f+Af  e  'l,t+4f  TV  e  /V  1  )i  2,t+kt  ~  \J2,t+At 

(2)  if  <-T  and  > F.(-T),  then 

ftSL  +  (1  -  <T*“  )F,  (,-T) 

It  can  be  shown  from  these  conditions  that  the  solution  is  plastic  when  rjTAt  — >  oo 


and  elastic  when  rjTAt  — >  0. 
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FIG.  2-12  Treatment  of  tension  cutoff 

2.3.2  THE  DUVANT-LIONS  TYPE  VISCOPLASTIC  CAP  MODEL 

The  viscoplastic  strain  rate  vector  and  hardening  parameter  are  respectively 
defined  as: 

&  =  -C -'[ct-ct] . (2.45) 

T 

$=-[k-k] . (2.46) 

T 

where  r  =  a  material  constant  called  the  relaxation  time;  the  pair  ( a,  k  )  =  the  stress  and 

hardening  parameter  of  the  inviscid  material  (a  bar  is  used  to  denote  the  variable  of  the 
inviscid  plastic  model)  which  can  be  viewed  as  a  projection  of  the  current  stress  on  the 
current  yield  surface;  k  and  $  =  hardening  parameter  and  its  differential  with  respect  to 
time. 


It  can  be  seen  from  Eq.  (2.45)  that  the  viscoplastic  strain  rate  is  simply  defined  by 
the  difference  between  the  true  stresses  and  the  stresses  obtained  by  the  inviscid  model 
which  is  quite  different  from  that  of  the  Perzyna  type  (Eq.  2.22). 
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2.3.2.1  STATIC  YIELD  FUNCTIONS 

The  Duvant-Lions  type  cap  model  plastic  yield  surface  function/  is  the  same  with 
the  Perzyna’s  type. 

2.3.2.2  SOLUTION  ALGORITHMS 

The  strain  rate  Eq.  (2.20)  and  (2.21)  are  integrated  over  a  time  step,  At,  from  t  to 
t+At,  to  yield  the  incremental  strains  and  stresses: 


As  =  As1'  +  A  svp . (2.47) 

A<r  =  CAse  =  c{As-As'v) . (2.48) 


where  As  =  the  total  incremental  strain  vector; 

A se  =  the  elastic  viscoplastic  incremental  strain  vector; 

A s'p  =  the  viscoplastic  incremental  strain  vector; 

Act  =  the  incremental  stress  vector. 

Based  on  the  Euler  method,  the  viscoplastic  incremental  strain  vector  A svp  can  be 
approximated  as 

. '(2.49) 

in  which  %  is  an  adjustable  integration  parameter,  0  <  y  <  1.  The  integration  scheme  is 
explicit  if  %  =  0  and  fully  implicit  if  %  =  1.  This  solution  algorithm  is  conditionally  stable 
when  x —  0-5  and  unconditionally  stable  when  />  0.5.  The  fully  implicit  integration 


scheme,  x=  1,  is  used  here  in  the  numerical  algorithm  just  for  simplification. 


Integrating  Eq.  (2.45)  over  a  time  step  At  gives 
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U  V>.,  a, .  | . (2.50) 

T 

Substitution  of  Eq.  (2.50)  into  Eq.  (2.48)  yields 

Aa  =  an+1-(Jn  =  C:Af -  — [crn+1-CT„+1] . (2.51) 

r 

By  solving  Ao;1+i  from  Eq.  (2.51),  one  obtains 
{a„  +  C:Ae)+  —  an+l 

cr  ,  = - A . . (2.52) 

"+1  ,  At 

1  +  — 

T 

where,  (cr„  +CAs)  may  be  treated  as  an  elastic  trial  stresses. 

Similarly,  we  obtain  the  hardening  parameter  may  be  expressed  as 
.  ^A t_- 

Kn  +  Kn+l 

K+l  = - . (2-53) 

1  +  — 

T 


The  numerical  algorithm  for  the  Duvant-Lions  viscoplastic  model  is  presented  in 
Table  2-2.  It  is  apparent  that  the  Duvant-Lions’  model  is  very  easy  implement,  since  the 
viscoplastic  solution  is  simply  an  update  of  the  inviscid  solution.  The  ease  of  numerical 
implementation  of  the  Duvant-Lions  model  is  apparent  compared  with  the  Perzyna  model, 
which  requires  many  matrix  operations. 


Table  2-  2  Numerical  algorithm  for  the  Duvant-Lions’s  viscoplastic  model 
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DATA  INPUT 

;  an,kn,Rs 

Trial  stresses 

:  =*„+CAe  ,kn 

If  f{<: ?.*.)<  o. 

elastic  crm+l  =  O'  ,  kn+l  =  kn  RETURN 

If  ,k„)<  0,  viscoplastic 

(a)  calculate  the  inviscid  solution:  ( <TIJ+1  ,£n+i ) 

(b)  update  to  viscoplastic  stress  and  hardening  parameter: 

At  _ 


(an  +C  :  Af)+  cr„ 

T 


k  +  —  k 


n+ 1 


1  + 


At 


k  a  1 1 


1  + 


At 


RETURN 


OUTPUT 


^/h-i  >  kn+i ,  £n+ 


2.4  ILLUSTRATION  EXAMPLE 

The  simulated  uniaxial  strain  test,  presented  by  Kantona  (1984),  was  used  to 
prove  the  adequacy  of  this  viscoplastic  cap  model  under  different  loading/unloading 
strain  rates. 

A  hypothetical  uniaxial  strain  loading  history:  the  axial  strain  of  the  soil  under 
compression  is  increased  at  a  constant  rate  (^,  =0.03  %/s)  for  1  second,  held  constant 
=0.0)  for  4  seconds,  unloaded  at  a  constant  rate  (^  =-0. 015%/s)  for  0.5  second,  and  held 
constant  afterwards  is  shown  in  FIG.  2-20. 

The  material  parameters  used  for  cap  model  are  those  for  McCormick  Ranch  sand 
given  by  Sandler  and  Rubin  (1979):  K  =  66.7  ksi;  G  =  40  ksi;  a  =  0.25  ksi;  {5  =  0.67  ksi’1; 
y=  0.18  ksi;  0=  0.0;  W=  0.066;  D  =  0.67  ksi’1;  R  =  2.5;  X0  =  0.189  ksi;  and  T=  0.0  ksi. 
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For  the  Perzyna  model,  the  two  parameters,  N  and/o,  were  assumed  to  be  1.0  and 
0.25ksi  based  on  experience  data,  respectively.  Three  values  of  the  fluidity  parameter 
(77  =  0.0035,  0.015  and  0.032)  were  examined  similarly.  According  to  Eq.  (2.22),  when  77 
decreases,  the  viscoplastic  strain  decreases,  and  the  stress  is  close  to  elastic,  which 
implies  the  axial  stress  will  increase.  The  stress  response  becomes  purely  elastic  as  77  — >  0, 
and  purely  plastic  as  77  — >  00. 

For  the  Duvant-Lions  model,  three  values  of  the  relaxation  time  (t  =  1.0,  0.25, 
0.125)  were  examined  to  illustrate  its  effects  on  the  stress  response.  As  shown  in  FIG.  2- 
21,  the  stress  response  increases  as  the  relaxation  time  rincreases.  According  to  Eq.  (2.45) 
and  (2.50),  when  rincreases,  the  viscoplastic  strain  decreases,  and  the  axial  stress  is  close 
to  elastic,  which  implies  the  stress  response  will  increase.  Although  it  is  not  plotted  in 
FIG.  2-21,  the  stress  responses  will  become  purely  elastic  as  r— >  00,  and  purely  plastic  as 
r— >  0. 


FIG.  2-13  Axial  strain  history  for  uniaxial  strain  test 
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FIG.  2-  14  Axial  stresses  for  different  values  of  r  and  rj 


By  comparing  the  stresses  resulting  from  the  two  models  in  FIG.  2-21,  it  can  be 
seen  that  each  pair  of  the  relaxation  time  and  fluidity  parameter  yields  nearly  the  same 
stresses.  For  instance,  the  axial  stress  history  with  r=  1.00  from  using  the  Duvant- Lions 
model  was  very  close  to  that  with  77  =  0.0035  from  using  the  Perzyna  model.  Likewise, 
stresses  obtained  from  using  the  Duvant-Lions  model  with  r  =  0.25  and  0.125  are  nearly 
the  same  as  those  obtained  from  using  the  Perzyna  model  with  77  =  0.015  and  0.032, 
respectively.  The  ratio  of  the  three  relaxation  times  is  8:2:1,  while  that  of  the  fluidity 
parameters  is  approximately  1:2:9.  Therefore,  a  certain  relationship  between  rand  77  may 
exist  and  the  viscoplasticities  of  these  two  types  may  be  equivalent  for  this  example. 
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2.5  MODEL  VALIDATOIN  AGAINST  EXPERIMENTAL  DATA 

Jackson  et  al.  (1980)  conducted  a  series  of  static  and  dynamic  tests  on  a  clayey 
sand.  These  tests  provided  data  for  validation  of  the  viscoplastic  cap  model  and  the 
associated  solution  algorithms. 

The  first  step  was  to  calibrate  the  material  parameters  for  yield  functions  and  the 
elastic  moduli  using  the  static  test  data.  The  static  test  data  consisted  of  the  stress  and 
strain  results  from  a  uniaxial  strain  test  and  two  triaxial  compression  tests  conducted  at 
confining  pressure  of  2.07  MPa  and  4.14  MPa,  respectively.  The  material  parameters 
obtained  to  fit  the  test  data  were:  K  =  2500  MPa;  G  =  1500  MPa;  a  =  3.654  MPa;  [3  = 
0.003  MPa1;  y=  3.500  MPa;  9=  0.263;  W  =  0.109;  D  =  0.05  MPa1;  R  =  1.5;  A0  =  0.3 
MPa;  and  T  =  0.0  MPa.  The  agreement  was  considered  to  be  good  both  qualitatively  and 
quantitatively. 

The  second  step  was  to  simulate  the  dynamic  stress-strain  relationship.  The  test 
data  were  obtained  from  dynamic  uniaxial  strain  tests,  each  of  which  was  conducted  at 
varying  strain  rate.  The  strain-histories  were  obtained  by  choosing  strain  and  time  values 
from  plots  of  vertical  stress  versus  time,  and  vertical  stress  versus  vertical  strain 
(Schreyer  and  Bean  1985).  The  maximum  strain  rate  in  the  dynamic  test  was 
approximately  200/s. 

The  additional  viscous  parameters  for  the  Perzyna’s  model  were:  ij  =  0.002msec"1; 


N  =  1.5;/o  =  1.0  MPa. 
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From  the  simulation  results  it  is  apparent  that:  (1)  formulations  are  capable  of 


predicting  the  dynamic  soil  response  well;  (2)  the  soil  responses  are  close.  For  the 
Perzyna  type,  a  normal  constant  strain  rate  of  0.0008/s  for  static  tests  was  used.  It  is 
apparent  that  the  soil  behavior  under  high  strain  rate  are  very  different  from  those 
obtained  in  static  test.  The  confined  modulus  and  the  strength  are  largely  increased  under 
high  strain  rate  loading.  The  viscoplastic  cap  models  capture  the  strain-rate  effects  very 
well. 


However,  there  are  some  slight  differences  between  predictions  of  the  two  models. 
For  instance,  the  initial  soil  stiffnesses  under  high  strain  rate  loading,  the  slopes  of  the 
responses  are  predicted  better  by  the  Perzyna’ s  model  than  those  by  the  Duvant-Lions’ 
model.  From  Eq.  (2.22)  and  (2.45),  the  Perzyna’s  viscoplastic  formulation  appears  to  be 
more  flexible  for  data  fitting  than  the  Duvant  -Lions’  formulation  due  to  more  viscous 
parameters  involved  (Tong,  X.,  and  Tuan,  C.Y.  2007).  Therefore,  the  Perzyna’s 
viscoplastic  cap  model  will  be  implemented  into  LS-DYNA  finite  element  code  to 
represent  the  soil  model  to  analyze  the  strain-rate  effect  due  to  explosion. 


CHAPTER  THREE  EQUATION  OF  STATES 
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3.1  INTRODUCTION 

An  ideal  liquid  or  gas  is  a  continuous  medium  with  neither  shear  or  frictional 
forces  acting  between  its  particles.  Hence  the  stress  at  a  given  point  does  not  depend  on 
the  orientation  of  the  small  surface  upon  which  it  acts.  In  actual  liquids  and  gases, 
frictional  forces  do  act  between  their  particles.  Solid  bodies  differ  from  liquids  and  gases 
in  that  they  transfer  shear  forces.  When  the  pressure  exceeds  a  certain  magnitude,  the 
bonds  between  the  particles  are  broken  so  the  material  is  compressed  and  the  solid  begins 
to  behave  like  a  fluid.  This  phase  change  depends  only  upon  the  magnitude  of  the 
pressure  and  the  temperature  (Grujicic  et  al.  2008).  The  state  of  a  medium  is  generally 
defined  by  a  combination  of  pressure  p  ,  density  p  ,  volume  V ,  temperature  T  ,  entropy  S, 
and  internal  energy  E.  All  these  quantities  are  related  by  thermodynamic  relations,  and 
only  two  of  these  quantities  are  independent.  The  general  form  of  P  =  P  (p  ,  E)  is  used 
herein  to  define  the  state  of  each  of  the  three  phases  of  the  soil. 

3.2  DEVELOPMENT  OF  SOIL  EQUATION  OF  STATES 

Any  equation  that  relates  the  pressure,  temperature,  and  specific  volume  of  a 
substance  is  called  an  equation  of  state.  There  are  several  equations  of  state,  some  simple 
and  others  very  complex.  Originally,  equation  of  states  were  mainly  used  in  physics  and 
thermodynamics,  an  equation  of  state  is  a  relation  between  state  variables.  More 
specifically,  an  equation  of  state  is  a  thermodynamic  equation  describing  the  state  of 
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matter  under  a  given  set  of  physical  conditions.  It  is  a  constitutive  equation  which 
provides  a  mathematical  relationship  between  two  or  more  state  functions  associated  with 
the  matter,  such  as  its  temperature,  pressure,  volume,  or  internal  energy.  Gradually, 
equations  of  state  are  found  that  are  useful  in  describing  the  properties  of  fluids,  mixtures 
of  fluids,  solids,  and  even  the  interior  of  stars. 

During  the  modeling  of  blast  loading  on  a  target  or  other  calculations  that  bring 
materials  together  at  high  velocities,  computer  simulations  of  materials  being  shocked  to 
high  pressure  and  then  releasing  to  low  pressure  are  performed.  Depending  on  the 
circumstances,  the  release  to  low  pressure  is  often  accompanied  by  release  to  a  very  low 
density.  Numerical  problems  leading  to  very  large  sound  speeds  or  to  negative  lagrangian 
volumes  have  been  encountered  during  numerical  simulation.  These  problems  can  be 
traced  to  the  behavior  of  the  equation  of  state  in  the  limit  as  the  density  becomes  much 
less  than  the  normal  or  reference  density. 

Since  all  three  phases  of  soil,  solids,  water  and  air,  have  significant  volume 
change  that  lead  to  change  pressure  and  density  under  blast  loading,  equations  of  state  are 
considered.  In  this  thesis,  which  is  focused  on  a  limited  number  of  equation  of  states  that 
can  be  used  for  solid  soil  finite  elements.  These  equations  of  states  include  Mie- 
Gruneison  equation  of  state,  Tillotson  equation  of  state  and  Kandaur  conceptual  equation  of 


state. 
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3.2.1  MIE-GRUNEISON  EQUATION  OF  STATE 

The  Mie-Gruneisen  equation  of  state  is  a  relation  between  the  pressure  and  the 
volume  of  a  solid  at  a  given  temperature.  It  is  often  used  to  determine  the  pressure  in  a 
shock-compressed  solid. 

If  the  pressure,  in  terms  of  energy  e  and  volume  v  is  expressed  as, 

P  =  f(e,v ) . (3.1) 


then  a  change  in  pressure  dPc an  be  written  as, 


dP  = 


(  dP) 

dv  + 

— 

UvJ 

e 

[de) 

de . 


.(3.2) 


Integration  of  this  equation  allows  the  pressure  to  be  expressed  in  terms  of  the  volume  v 
and  energy  e  relative  to  the  pressure  at  a  reference  volume  v0and  reference  energy  eQ. 


The  integration  can  be  performed  along  any  path  desired  and  it  is  convenient  to 
integrate  first  at  constant  energy  from  v0  to  v  ,  and  then  at  constant  volume  from  e0  to  e  , 
giving, 


P  = 


dP 

de 


de 


(3.4) 


The  Gruneisen  Gamma  is  defined  as, 
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(3.5) 


and  if  it  is  assumed  that  r  is  a  function  of  volume  (or  density),  only  then  can  the  second 
integral  above  be  evaluated, 


The  first  integral  is  a  function  only  of  volume  and  the  reference  energy  q, .  If  the 
reference  state  is  denoted  by  er ,  then  since, 


j(ff)  'M.l  V . (3.7) 

The  equation  becomes, 

P  =  Pr(v)+IH[e-er{v)\ . (3.8) 

v 


This  equation  is  generally  known  as  the  Mie-Gmneisen  form  of  equation  of  state.  In  LS- 
DYNA,  it  can  be  expressed  as, 

PoC'fi  1+  — 

p  =  r - - - - - +  (r0  +  aju)E . (3.9) 

1-(5,-1)//-52^--53-^t 

L  A  +  1  (//  + 1)  J 

Where  C  is  the  intercept  of  the  Shock  velocity-Particle  velocity  curve;  Sr  S2  and  S 3  are 
the  coefficients  of  the  slope  of  the  Shock  velocity-Particle  velocity  curve;  y()  is  the 

Gruneisen  gamma;  a  is  the  first  order  volume  correction  to  y{) ;  and  //  =  --!. 
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3.2.2  TILLOTSON  EQUATION  OF  STATE 

This  form  of  equation  of  state  (Tillotson  1962)  was  derived  to  provide  a 
description  of  the  material  behavior  of  solid  elements  over  the  very  wide  range  of 
pressure  and  density  encountered  in  hypervelocity  phenomena. 

Not  only  must  such  an  equation  of  state  describe  normal  density  material  and  its 
condition  after  shock,  but  also  its  expansion  and  change  of  phase  in  cases  where  the 
shock  energy  has  been  sufficient  to  melt  or  vaporize  the  material.  The  pressure  range  can 
be  so  large  that  the  “low  pressure”  regime  of  this  form  of  equation  of  state  is  defined  as 
from  0  to  10  Mbar  and  “high  pressure”  from  10  to  about  1000  Mbar.  Thus  any  pressure 
and  results  from  normal  laboratory  experiments  cover  only  the  “low  pressure”  region.  For 
the  derivation  of  an  equation  of  state  for  the  “high  pressure”  region,  analytic  forms 
provide  best  fit  interpolations  between  Thomas-Fermi-Dirac  data  at  high  pressures  (above 
50  Mbar)  and  experimental  data  at  low  pressures.  The  formulation  is  claimed  to  be 
accurate  to  within  5%  of  the  Hugoniot  pressure  and  to  within  10%  of  the  isentrope 
pressures.  It  is  therefore  a  very  useful  form  of  equation  of  state  for  hypervelocity  impact 
problems. 

The  total  range  of  the  pressure-volume  plane  is  divided  into  four  regions  as  shown 


in  the  FIG.  3-1. 
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FIG.  3-  1  Regions  of  Interest  in  the  (p,  v)  Plane 

The  region  to  the  left  of  the  Hugoniot  can  only  be  reached  by  adiabatic  (non¬ 
shock)  compression  and  is  not  relevant  to  impact  problems.  It  is  therefore  excluded  from 
the  formulation.  Region  I  represents  the  compressed  phase  of  the  material  and  extends 
vertically  to  shock  pressures  of  about  1000  Mbar.  Region  II  describes  material  which  has 
been  shocked  to  energy  less  than  the  sublimation  energy  and  will  therefore,  on  adiabatic 
release,  returns  to  zero  pressure  as  a  solid.  There  is  no  provision  in  this  equation  of  state 
to  describe  the  material  at  pressures  less  than  zero.  Region  IV  is  the  expansion  phase  of 
material  which  has  been  shocked  to  energy  sufficiently  large  to  ensure  that  it  will  expand 
as  a  gas.  For  large  specific  volumes,  the  formulation  for  Region  IV  extrapolates  to  an 
ideal  gas  limit.  It  is  desirable  or  even  necessary,  to  ensure  that  the  formulations  in  each 
region  provide  continuous  values  of  the  pressure  and  its  first  derivatives  at  the  boundaries 
between  contiguous  regions.  Region  III  lies  between  Regions  II  and  IV.  In  this  region 
the  pressure  is  calculated  as  a  mean  between  that  calculated  for  Regions  II  and  IV. 


Let: 
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77  =  — . (3.10) 

Po 

M  =  rj~  1 . (3-11) 

coQ=l  +  -^ . (3.12) 

eoB 


where  p  is  the  density,  pQ  is  the  reference  density,  e  is  the  energy  and  e0  is  the 
reference  energy. 

For  Region  I  (p  >  0)  the  pressure  Pl  is  given  by  a  Mie-Gruneisen  equation  of 
state  but  since  the  formulation  is  to  be  valid  for  a  very  large  range  of  pressure,  the 
Gruneisen  Gamma  is  a  function  of  both  v  and  e ,  not  just  a  function  of  v  alone.  The 
constants  fit  the  low  pressure  shock  data  but  they  are  adjusted  to  fit  the  asymptotic 
Thomas-Fermi  behavior  for  the  variation  of  pressure  at  maximum  compressions  (like  a 
monatomic  gas).  The  formulation  for  Region  II  is  as  for  Region  I  with  a  slight 
modification  to  one  term  to  avoid  difficulties  as  m  goes  increasingly  negative.  In  Region 
IV  the  formulation  is  chosen  to  give  the  correct  behavior  both  at  high  pressure/normal 
density  and  for  very  large  expansion  (where  it  must  converge  to  an  ideal  gas  behavior). 
With  these  constraints  the  different  formulations  are  given.  For  region  I  (p  >  0)  ,  the 
pressure  P1  is  given  by, 


P>  = 


CL  +  - 


CO, 


rjp0e  +  Ap  +  Bp2 


0  ) 


V 


(3.13) 


For  region  II  (//  <  0,  e  <  es )  ,  the  pressure  P2  is  given  by, 
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P2  = 


'  b ' 

CL  H - 

V  Mq  J 


i!pQe  +  Afu . (3.14) 


For  region  III  (//  <  0,  es  <  e  <  es )  ,  the  pressure  P3  is  given  by, 
(P,-P2\e-e,) 


P  =  P  + 
1  3  1  2  ^ 


((3  <0 


.(3.15) 


For  region  IV  (//  <  0,  e  >  es )  ^  the  pressure  P,  is  given  by, 


P4  =  ar/p0e  + 


brjp^e 


V 


+  Ape 


P* 


1 


.(3.16) 


where  x  =  1 - .  In  the  Tillotson  equation  of  state,  a,  b  ,  A,  B  ,  e0,  e9,  and  es  ai'e 

V 


constants. 


The  Mie-Gruneisen  equation  of  state  and  Tillotson  equation  of  state  can  be  used 
for  soil  behavior  simulation  model  and  ensure  a  unique  solution.  However,  the  limitation 
is  soil  with  Mie-Gruneisen  equation  of  state  or  Tillotson  equation  of  state  is  defined  as  a 
unit  material  and  leads  to  a  simplified  bulk  modulus  and  mechanical  pressure  in  the 
calculation  process. 


3.2.3  MURRAY’S  EQUATION  OF  STATE  FOR  UNSATURATED  SOILS 

The  prediction  of  soil  behavior  is  intrinsically  linked  to  the  need  to  determine  the 
controlling  stresses  in  the  soil.  For  saturated  soils,  Terzaghi  (1936)  proposed  an  equation 
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for  effective  stress  which  controls  the  shear  resistance  and  volume  changes.  The  effective 
stress  can  be  written  as 

p'  =  P~uw . (3.17) 

where  p'  is  Terzaghi’s  mean  effective  stress,  p  is  the  mean  total  stress  and  uw  is  the 
pore-water  pressure. 

The  concept  of  the  stress  state  variable  (p-uw)  controlling  the  behavior  of 
saturated  soils  has  proven  very  useful  and  has  been  shown  to  be  valid  in  practice.  For 
unsaturated  soils,  however,  the  search  for  a  reliable  stress  state  variable  equation, 
independent  of  soil  properties,  has  proven  unsuccessful.  As  described  by  Fredlund  and 
Rahardjo  (1993),  a  number  of  such  equations  have  been  proposed.  The  original 
suggestion  of  Bishop  (1959)  can  be  written  as 

p'b  ={p-ua)+  x{ua  -«w) . (3.18) 

where  p'B  is  Bishop’s  mean  effective  stress,  ua  is  the  pore-air  pressure  and  %  is  an 
empirical  parameter. 

A  major  obstacle  to  the  use  of  Eq.  (3.18)  lies  with  the  parameter  %  .  This  is  usually 
ascribed  the  range  of  values  0  <  %  <  1  and  has  been  shown  to  be  dependent  on  the  stress 
path  and  the  process  to  which  the  soil  is  subjected  (Jennings  and  Burland  1962;  Blight 
1965;  Morgenstem  1979). 

Although  it  is  desirable  that  the  concept  of  effective  stress  for  saturated  soils 
extended  to  unsaturated  soils  and  that  soil  properties  such  as  the  volumes  of  the  various 
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phase  (solid  particles,  water  and  air)  are  not  included  in  any  formulation  of  controlling 


stresses,  experiments  have  demonstrated  the  inadequacy  of  any  such  relationship.  For  this 
reason,  researchers  have  turned  to  examining  the  use  of  the  independent  stress  state 
variables  (p-ua),  (p-uw)  and  (ua  -uw) to  describe  the  mechanical  behavior  of  soils. 

Fredlund  and  Morgenstern  (1977)  concluded  from  theoretical  considerations  that  any  two 
of  these  three  stress  state  variables  can  be  used  to  describe  the  behavior  of  an  unsaturated 
soil.  However,  there  are  inconsistencies  in  experimental  results  not  readily  answered  by 
constitutive  modeling  using  these  parameters  (Wheeler  and  Sivakumar  1995).  A  logical 
interpretation  of  experimental  data  is  essential  to  an  appreciation  of  soil  behavior,  and  a 
clear  pricture  does  not  always  emerge  using  independent  stress  state  parameters,  as  these 
interact  in  response  to  external  stimuli.  In  this  respect,  it  appeal's  that  the  volumes  of  the 
phases  play  an  important  role  in  controlling  the  stresses  in  unsaturated  soils,  and  this  is 
demonstrated  in  the  analysis  and  the  comparisons  with  both  consolidation  data  and 
critical  state  data  which  follow. 

Murray  (2002)  examined  the  significance  of  the  relative  volumes  of  the  phases, 
and  the  interactions  between  the  phases,  on  the  stress  regime  under  equilibrium 
conditions.  First,  a  description  of  the  significance  of  enthalpy  in  soils  relating  pressures, 
volumes,  and  internal  energy  sources  is  presented,  followed  by  an  examination  of 
Terzaghi’s  effective  stress  equation  in  terms  of  the  enthalpy  of  a  saturated  system.  This 
approach  is  then  extended  to  unsaturated  soils  to  develop  an  equation  of  state  that 
includes  the  average  volumetric  “coupling”  stress  p'c  .  This  links  the  stress  state  variables 


and  the  volumes  of  the  phases. 
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The  general  equation  of  state  for  unsaturated  soil  can  be  expressed: 


P  =  uwnw  +  b(nw  +  ns ){ua  -uw)+uana  +auwns  +  p[. . (3.19) 

P'c  =  {p  ~  u a )  +  sinw  +  ns ) . (3-20) 


where  na  (na  =  Va  /V)  ,  nw  (nw  =VW/V)  and  ns  (ns  =VS/V)  are  the  volume  fraction  of 

air,  water  and  solid  phase  respectively,  a  is  a  dimensionless  parameter  with  a  minimum 
value  of  1,  b  is  a  dimensionless  number  influenced  by  the  structure  and  size  of  the 
saturated  packets  and  5  is  the  suction.  (nw+ns)  represents  the  total  volume  of  the 
saturated  packets  per  unit  volume  of  soil.  Using  Eq.  (3.19)  it  is  possible  to  highlight  the 
significance  of  the  stress  state  variables  (p-ua)  ,  (p-uw)  ,  and  (ua  —  uw )  for 
unsaturated  soils  and  their  implicit  relationship  with  the  volumes  of  the  phases. 

FIG.  3-2  and  FIG.  3-3  have  been  prepared  based  on  the  experimental  data 
reported  by  Wheeler  and  Sivakumar  (1995)  and  Toll  (1990). 


FIG.  3-  2  Variation  of  specific  volume  during  ramped  consolidation  at  different  suction 

(Murray,  2002) 
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FIG.  3-3  q  versus  p'c  at  critical  states 


(MuiTay,  2002) 


However,  the  Murray’s  equation  of  state  has  the  limitations  and  shortcomings:  (1) 
there  are  no  enough  experimental  data  at  present  to  define  precisely  the  suction- 
dependence  of  material  parameters  and  (2)  the  average  volumetric  coupling  stress  p'c 

represents  the  microscopic  forces  between  particles.  Under  high  strain  rate  loading,  like 
blasting  loading,  the  average  volumetric  coupling  stress  doesn’t  play  an  important  role. 


3.3  KANDAUR’S  CONCEPTUAL  MODEL  OF  EOS 

Soils  are  composed  of  particles  of  various  materials-  called  phases.  The  majority 
of  the  solid  mineral  particles  consists  of  silicon  which  can,  therefore,  be  taken  as 
representative,  the  other  water  and  air. 
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Let  As ,  Aw  and  Aa  denote  the  relative  volume  of  the  solid  particles,  water  and  air, 


respectively,  i.e.  the  volume  of  the  corresponding  phase  in  a  unit  volume  of  soil;  then 

As+Aw+Afl=l  . (3.21) 

The  quantities  ps ,  pw  and  pa  are  the  material  densities  of  each  phase  and  p{)  is  the 
initial  density  of  the  soil  as  a  whole.  We  then  have 

A)  =  APs  +  AwPw  +  A,  Pa  . (3-22) 

In  soils,  two  deformation  mechanisms  exist: 

a)  at  low  pressures,  the  soil  skeleton  deformation  is  determined  by  the  elastic 
deformations  of  bonds  on  the  contact  surfaces  of  grains  and,  at  high  pressures,  it 
is  determined  by  a  failure  in  bond  and  displacements  of  the  grains  (plastic 
deformation); 

b)  the  deformation  of  all  the  soil  phases,  determined  by  their  volume  compression. 
When  the  soil  is  being  compressed,  both  mechanisms  are  always  acting 
simultaneously.  At  certain  phases  of  the  loading  process,  however,  one  of  the 
mechanisms  predominates  to  such  a  degree  that  the  other  may  be  neglected. 

A  dry  soil  contains  air  and  a  small  amount  of  water,  whose  compressibility 
considerably  exceeds  that  of  the  skeleton;  therefore,  with  static  and  dynamic  loading,  the 
first  mechanism  becomes  influencial  while  the  other  is  negligible;  with  increasing 
pressure,  the  gain  bonds  are  deformed  and  displaced  and  the  soil  is  compacted  so  that  the 
second  mechanism  becomes  more  and  more  effective  until  it  reaches  a  definite 
overbalance,  while  the  first  becomes  negligible.  The  dependence  of  pressure  on  the 
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relative  volume  deformation  is,  for  this  case,  plotted  in  FIG.  3-4  (Henrych  1979).  The 
second  mechanism  predominates  for  a  >  <rB. 


Relative  volume  deformation  0 

FIG.  3-  4  Relationship  between  stresses  and  relative  volume  deformation  for  solids 

In  a  saturated  soil  the  salts  on  the  grain  contacts  are  dissolved  and  the  bonds 
weakened.  Under  a  rapid  dynamic  loading,  the  water  and  air  have  a  higher  resistance  than 
the  contact  bonds  of  the  skeleton  grain.  The  deformation  and  resistance  of  the  soil  are 
determined  by  the  dominating  second  mechanism,  particularly  by  the  water  and  air 
deformation;  the  solid  phase  becomes  effective  only  at  high  pressures  (hundreds  and 
thousands  of  kp/cm2).  The  relationship  between  cr(0)  and  volume  deformation  under  this 
situation  is  shown  in  FIG.  3-5.  However  under  a  slow  static  loading  of  the  saturated  soil, 
the  water  and  air  are  pressed  out  of  the  void  and  the  compressibility  is  mainly  given  by 


the  solid  skeleton  compressibility. 


o 

ysy_ 

Relative  volume  defomiation  0 

FIG.  3-  5  Relationship  between  stresses  and  relative  volume  defomiation  for  liquids, 

gases,  etc 
(Henrych,  1979) 

The  diagram  of  a  block  grain  medium,  according  to  I.I.  Kandaur  (Henrych,  1979), 
is  illustrated  in  FIG.  3-6.  The  cavities  between  blocks  are  filled  with  water  and  air. 
Between  the  comers  are  elastobrittle  bonds.  With  loading,  the  medium  deformations 
consist  of  the  deformations  of  the  elastobrittle  bonds,  which  are  disturbed  with  a 
simultaneous  mutual  displacement  of  the  static  blocks  (first  mechanism)  and  the  void 
filled  with  water  and  air  (second  mechanism).  The  forces  of  the  elastobrittle  bonds  and 
the  forces  of  friction  between  the  solid  particles  act  within  the  scope  of  the  first 
mechanism.  The  forces  depend  on  the  volume  change  of  the  individual  phases  then  act 
within  the  range  of  the  second  mechanism.  With  fast  dynamic  deformation  the  water  and 
air  are  cannot  escape  from  the  cavities  through  the  spaces  between  the  blocks;  with  a 
slow  static  defomiation  the  water  and  air  are  forced  through  the  spaces  between  the 
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blocks  into  less  loaded  surroundings  and  the  dominant  resistance  is  offered  by  the  bonds 


between  blocks  and  by  the  blocks  themselves. 


FIG.  3-  6  Schematic  representation  of  a  block  grained  medium  with  elastobrittle 

linkages  between  the  blocks 
(Henrych,  1979) 


The  medium  shown  in  FIG.  3-6  corresponds  to  the  rheological  model  illustrated 
in  FIG.  3-7,  which  covers  both  mechanisms  and  applies  to  a  dynamic  loading  (water  and 
air  are  not  forced  out  of  the  voids).  This  model  is  used  to  derive  the  equation  of  state  for 
the  adiabatic  process.  With  small  pressure  and  dry  soils  the  first  deformation  mechanism 
is  a  decisive  factor  as  it  corresponds  to  the  elements  D,  E,  i.e.  to  the  grain  friction 
proportional  to  normal  pressure,  and  to  the  resistance  of  the  crystal  bonds,  which  is 
represented  by  a  series  of  filaments  stretched  and  broken  as  the  deformation  develops. 
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FIG.  3-  7  Schematic  diagram  of  a  rheological  model  of  the  medium 

(Henrych,  1979) 


With  water-bearing  soils  and  for  higher  pressures  with  dry  soils,  the  second 
deformation  mechanism  represented  by  the  elements  A,  B,  C  predominates.  Obviously, 


P  =  Pa+Pb+Pc  . (3.23) 

V  =FS  +  Fp  . (3.24) 

Vp=Vw+Va  . (3.25) 

V.=A,V0  . (3.26) 

VW=AWV0  . (3.27) 

Va=\V0  . (3.28) 


where  Pa  ,  Ph  and  Pc  are  the  forces  in  branches  a,  b  and  c,  respectively;  V  is  the  soil 
volume,  V0  is  the  soil  initial  volume  and  V  is  the  void  volume. 


From  equations  (3.23)  to  (3.28),  we  obtain 


dv  =  dVs  +  dVn 

A  y 
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.(3.29) 


dP  =  —dVn  +  —  +^dVP 


.(3.30) 


dVP  = — —dPh  + — -dPh 

dPh  dPu 


.(3.31) 


and  hence 


(  dV  )  [  dV  dV  )  dP  dP 
dP-  dV - -dP  — ^  +  — a-  +  =0 


8Pb  dPh)  8Vp  dVp 


.(3.32) 


Then  dependence  of  the  loading  on  deformation  in  phases  1  and  2  is  given  by  the 


Hooke  law,  so  that 


^  A,V0 


.(3.33) 


AwV0 


.(3.34) 


where  k*  ,  k*  are  the  coefficients  of  volume  deformation  of  the  mineral  skeleton 


particles  and  of  water,  respectively. 


In  element  C  holds  the  equation  of  state  of  a  polytropic  gas,  which  can  be  written 
in  the  form 


''  '•  '<J'>  . (3-35) 


where  P0  is  the  atmospheric  pressure,  aa  is  a  constant  and  k  is  the  coefficient  of 


adiabaticity.  Then, 
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(3.36) 


In  the  element  a,  the  relationship  between  loading  and  deformation  is  determined 
by  the  dry  friction  produced  by  a  force  P’  between  the  blocks,  proportional  to 


deformation: 

Pa=fP'  . (3-37) 

P'  =  KPAVP  . (3.38) 

A Vp=VP-(Aw+Aay0  . (3.39) 

where  /  is  the  coefficient  of  friction  of  the  mineral  particles  and  Kp  is  the  coefficient  of 
proportionality.  From  equations  (3.37)  to  (3.39),  follows  the  coefficients, 

Pa=<P*VP  . (3.40) 

V  =  Kpf  . (3.41) 

which  are  constant  for  a  given  soil  and  moisture,  so  that, 


The  force  in  each  filament  of  the  element  E  obeys  the  Hooke  law  until  the 
filament  breaks.  But  the  strength  of  the  individual  filaments  is  different  and,  therefore, 
the  force  Pc  in  the  aim  c  is  expressed  as, 


Pc  =  EAVP  . (3.43) 

where  E  is  a  variable  deformation  modulus,  which  may  be  written, 

E  =  E0(\-AE)  . (3.44) 


where  E0  is  a  constant.  With  regard  to  the  statistical  law  of  disturbance, 
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AE  =  \e~*dx  . (3.45) 

x  =  -BAVp  . (3.46) 

We  obtain 

Pc  =  EQAVpeBAVp  . (3.47) 

so  that, 

^  =  -E0(\  +  BAViyAV"  . (3.48) 


Substituting  the  equation  (3.48)  into  equation  (3.32),  the  equation  can  be  obtained, 


dp  + 

{dV-d4\ 

k  ks  y 

V  ) 

+  EQ(l  +  BAVP)eBAVp  -cp 


0  . (3.49) 


For  the  initial  condition  (3.50),  it  is  possible  to  obtain  the  solution  of  equation  (3.49)  in 
the  form  (3.51).  Because  of  the  their  inordinate  complexity,  neither  equation  (3.49)  nor 
its  solution  have  as  yet  been  used  for  dynamic  problems,  even  if  it  determines  the 
behavior  of  soil  with  sufficient  accuracy. 


W  )  . (3.50) 

P(v)  =  P  . (3.51) 


For  the  solution  of  soil  dynamics  problems  the  equation  of  state,  derived  by  G.M. 
Lyakhov  (Henrych,  1979),  is  more  suitable.  This  equation  is  based  on  the  second 
mechanism  of  soil  deformation,  i.e.  the  volume  compression  of  all  phases;  in  deriving  it 
Lyakhov  stalled  from  the  equations  of  state  of  the  individual  phases. 


For  air,  the  equation  of  state  can  be  expressed  in  the  form, 
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P.  =  Pa 


(  \k“ 

Pa 

\PaO 


.(3.52) 


where  P0  is  the  atmospheric  pressure,  it  can  be  expressed, 


P0  = 


PaOCaO 


.(3.53) 


pa0  is  the  density  of  air  at  atmospheric  pressure,  ca0  is  the  velocity  of  sound,  is  the 
density  of  air  at  pressure  and  ka  is  the  exponents  of  the  specific  entropy  of  the  air. 


For  water,  the  equation  of  state  can  be  expressed  in  the  form, 


Pw  =  Po  + 


/^wO^*wO 


f  p  ^ 
V  PwO  J 


-1 


.(3.54) 


For  solid,  the  equation  of  state  can  be  expressed  in  the  form, 


Ps=P0  + 


PsQCsO 


f  Ps  ^ 


\PsOj 


-1 


.(3.55) 


These  pai'ameters  of  equations  of  state  are  summarized  in  Table  3-1. 


Table  3-  1  Equation  of  state  pai’ameters  for  saturated  soil 


p0  (kg/m3) 

c^(km/s) 

k 

Air 

1.2(ao) 

0.34  (ca0) 

1.4  (ka) 

Water 

1000(pwo) 

1.50  (cw0) 

UK) 

Solid 

2650(ao) 

4.50  (cs0) 

3  (ks) 
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For  solid,  water  and  air,  the  relative  volume  by  Aa  ,  Aw  ,  As ,  the  density  by  pa0  , 
Pwo  ’  Pso  ’  and  the  velocity  of  sound  by  ca0  ,  cw0  ,  cs0  ,  respectively,  at  an  initial 
(atmospheric)  pressure  p=po-  Because  of  the  different  compressibilities  of  the 
components,  their  relative  content  at  pressure  p  will  be  different  from  that  at  pressure 
p-po-  If,  at  pressure  p,  the  content  of  the  components  is  denoted  by  A*  ,  A*  ,  A* ,  the 
specific  volume  by  Va  ,  V  ,  Vs  and  the  soil  density  by  p  ,  it  follows  from  equation  (3.52) 


that, 


rAy- 


\AaJ 


f  \ 


v. 


a  0 


J 


.(3.56) 


It  can  be  rewritten  as. 


A  =  A 

a  a 


r  P  ^  ~k ' 


V  Po  J 


.(3.57) 


Similarly,  for  water 


kJPK  -P»)  !  ) 


—k~ 


V  Pw 0  Cw 


.(3.58) 


For  solid  paiticles 


A,.  =  A 


AriA»)+1 


PsO  Cs 


.(3.59) 


Because  the  density  increments  have,  due  to  compressibility,  the  density  of  a 


three-phase  medium  at  pressure  p  will  be, 


p  -  A)  (A  +  ) 
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(3.60) 


Thus,  water-bearing  and  dry  soils  may,  within  a  certain  pressure  range,  be 
considered  as  three-phase  media.  The  smaller  the  value  of  Aa  and  the  greater  the  value  of 

Aw  in  the  soil  voids,  the  lower  the  pressure  Pmin  corresponding  to  the  lower  limit  of 
applicability  of  this  model.  For  water-bearing  soils  Pmw  =  PQ  when  Aa  =  0  and 
Pmin  =  500  to  800  kp  when  Aa  =  0.04  to  0.05.  For  dry  soils  with  Aa  =  0.3  to  0.4,  the 
value  ofPmin  increases  up  to  several  hundred  to  several  thousands  of  atmospheres.  The 
upper  limit  is  bounded  by  the  validity  limits  of  the  equations  of  state  of  the  individual 
components. 

3.4  USER  DEFINED  EQUATION  OF  STATE 

To  improve  simulation  results,  an  equation  of  state  was  defined  for  LS-DYNA 
dynamic  simulation  software. 

The  conservations  of  mass,  momentum  and  energy  in  a  soil  medium  from  the 
initial  state  (denoted  by  the  subscript  0)  to  the  state  under  shock  loading  (denoted  by 
subscript  H )  are  expressed  by  (3.61),  (3.62)  and  (3.63),  respectively: 


Po  us  =p(Us  ~up ) . (3.6i) 

PH  =  Pi)  us  UP  . (3-62) 

E„-E„=!f{Vo-V„) . (3.63) 


where  U s  is  the  shock  velocity,  and  up  is  the  particle  velocity. 
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A  series  of  plate  impact  experiments  were  performed  on  a  soil  at  various  levels  of 
water  saturation  by  Chapman  et  al.  (2006).  The  Hugoniot  was  determined  using  a 
reverberation  technique.  The  Hugoniot  is  presented  in  terms  of  the  measured  shock 
velocity  and  particle  velocity  in  FIG.  3-8,  and  in  terms  of  stress  and  particle  velocity  in 
FIG.  3-9.  The  densities,  degrees  of  saturation  and  shock  wave  velocities  in  the  soil 
specimens  are  summarized  in  Table  3-2. 


•  22%  sat  1.84  g  cm-3 

♦  20%  sat  1.81  g  cm-3 
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FIG.  3-  8  Shock- velocity  vs.  particle-velocity 
(Chapman  et  al.  2006) 
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FIG.  3-  9  Stress  vs.  particle- velocity 
(Chapman  et  al.  2006) 


Table  3-  2  Plate  impact  test  data 


Moisture,  % 

0 

10 

20 

22 

Saturation,  % 

0 

32 

64 

70 

Density,  kg  m-3 

1430 

1530 

1810 

1840 

Shock  velocity,  km/s 

1.44 

1.45 

1.90 

2.68 

Hugoniot  curves  are  often  expressed  as  a  relation  between  shock  velocity  and 
particle  velocity  by  least-square  curve  fitting  the  shock  loading  data  (Zukas  1990).  For 
many  materials,  the  Hugoniot  can  be  expressed  as  a  linear  relation  between  shock 


velocity  Us  and  particle  velocity  u ,, : 
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U s  =  C0  +  s  Up  . (3.64) 

where  C0  is  the  sound  speed  at  ambient  pressure  and  temperature,  and  5  is  the  slope  of 
the  lineal-  relation,  both  obtained  experimentally. 

Dividing  both  sides  of  (3.64)  by  U s  yields 

1  =  —  +  5— . (3.65) 

Us  Us 

From  (3.61),  the  volumetric  strain  can  be  expressed  as 

—  =  1-  ^  =  A  . (3.66) 

Us  P  Vo 

Substituting  (3.66)  into  (3.65)  yields 

r.  ,  (3.67) 

1-5  A 

From  (3.66), 

uP  =  US  a  =  _^°A. . (3.68) 

1-5  A 


Let 

P  =  —  -1  =  ^-1 . (3.69) 

Po  V 

Substitute  (3.66)  into  (3.69), 

A  =  . (3.70) 

1  +  jU 

Substituting  (3.67),  (3.68)  into  (3.62)  yields 


6  ) 

f  C0  A  ) 

I  _  Po  Co  A 

V 1  —  5  A  J 

V1  —  5  A  J 

<N 

<f 

1 

1 

PH  =  Pi) 


(3.71) 


Substituting  (3.70)  into  (3.71)  yields 
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A,  Q 


Pu  = 


p 

1  +  Pj 


\-s 


P 


(3.72) 


V  1  +  P) 

Plate  impact  experiments  have  been  conducted  by  many  researchers  to  provide 
Hugoniot  data  for  various  materials.  Jones  and  Gupta  (2000)  conducted  shock  wave 
experiments  to  determine  the  refractive  index  and  shock  velocity  of  quartz. 
Braithwaite  et  al.  (2006)  obtained  the  shock  Hugoniot  properties  of  quartz  feldspathic 
gneiss  by  plate  impact  experiments.  The  relationship  between  shock  velocity  U s  and 


particle  velocity  up  of  solid  can  be  obtained  from  FIG.  3-10. 


0.1  0.15  0.2  0.25  0.3  0.35 

Particle  Velocity  (km/s) 


FIG.  3-  10  Shock-velocity  dependence  on  particle-velocity  for  quartz 


Us  =  6.319  +  1.41nP . (3.73) 

Nagayama  et  al.  (2002)  obtained  a  linear  relation  between  the  shock  velocity  and 
particle  velocity  of  water  from  high  velocity  impact  tests,  as  presented  in  FIG.  3-11. 
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FIG.  3-11  Shock-velocity  dependence  on  particle-velocity  for  water 


Us  =1.460  +  2.0m  p 


(3.74) 


Kim  et  al.  (1991)  investigated  the  Hugoniot  data  of  dry  air  and  derived  an 
expression  for  adiabatic  exponent  for  shock  compressed  dry  air  in  FIG.  3-12. 


FIG.  3- 12 


Shock-velocity  dependence  on  particle- velocity  for  air 
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(3.75) 


Us  =0.241  +  1.06 uP 

These  Hugoniot  data  for  quartz  sand,  water  and  air  are  used  in  the  equations  of 
state  and  are  summarized  in  Table  3-3. 


Table  3-  3  Equation  of  state  parameters  for  soil 


A 

p  (kg/m3) 

C0  (km/s) 

s 

k 

To 

Solid 

0.7  (As) 

2650(a) 

6.319 

1.41 

3  (ks) 

1.0 

Water 

0.2(AW) 

1 000 ( Av) 

1.460 

2.00 

7  (kw) 

0.6 

Air 

0.1(Aa) 

1.2(A) 

0.241 

1.06 

1.4  (ka) 

0.0 

Dry  soil 

1.0 

1430 

0.530 

1.64 

— 

0.11 

Saturated  soil 

1.0 

1840 

0.320 

4.92 

— 

0.11 

An  equation  of  state  for  states  more  general  than  the  uniaxial  strain  condition  in 
the  plate  impact  experiments  can  be  expressed  as  (Zukas  1990): 

P  =  py(V)E  =  ^p-E . (3.76) 

where  /(’ V )  is  the  Gruneisen  parameter,  and  E  is  internal  energy  per  unit  mass.  If  shock 
pressure  PH  and  internal  energy  EH  are  associated  with  a  specific  volume  V  from  a 
Hugoniot  curve,  the  shock  pressure  is  expressed  as 


PH=py{V)EH  =^Eh . (3.77) 

If  the  Hugoniot  is  the  reference  state,  the  equation  of  state  can  be  expressed  as 
P~PH  =^p-(E-EH)  . (3.78) 


Substituting  (3.63)  into  (3.78)  yields  the  equation: 
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p  =  pk+^(e-e„)-~  (va-v)  =  p„[i-^y^(E-Ea). 


Let 


r  =  «  +  (/()  -  a) 


v_ 

vn 


or 


y  -  a  _  y0  -  a 


V 


Vn 


Substituting  (3.69)  into  (3.80)  yields 


Y  =  vr(r<>  +  an)  =  7^—  O'o  +  an) 


V 


i  +  n 


The  internal  energy  per  unit  initial  volume  is: 
E 


Ev  =■ 


V 


Substituting  (3.72),  (3.82)  and  (3.83)  into  (3.79)  yields 


A)  C0  // 


P  = 


1  + 

a  , 

V  2  j 

2 

+ 


(y0  +  a//)£K. 


.(3.79) 


.(3.80) 


.(3.81) 


.(3.82) 


.(3.83) 


.(3.84) 


(1  +  //-S//)2 

where  the  initial  internal  energy  E0  in  (3.79)  corresponds  to  the  mechanical  work  done 
by  the  hydrostatic  pressure  in  soil  due  to  gravity.  Using  the  parameters  given  in  Table  3- 
3,  Equation  (3.84)  can  be  used  to  calculate  the  pressures  in  the  three  phases  of  the  soil. 


For  solid,  the  equation  of  state  can  be  expressed  as, 


p  =  (2650X6.319)7  (1  +  0.5//)  +  , 

(1-0.41//)2 


(3.85) 


For  water,  the  equation  of  state  can  be  expressed  as, 
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P  =  (1  OOOXl  -460)2  //  (1  +  0.7//)  +  (0  6)E,. 

(1  -m)1 

For  air,  the  equation  of  state  can  be  expressed  as, 

„  (1-2)(0.241)2  jujl  +  ju) 

(1-0.06  pf 

The  bulk  modulus  of  the  soil  can  be  calculated  as, 


(3.86) 


(3.87) 


„  r 2 
Po  Co 


K  = 


1  + 


'  v  A 

l_Zo 

v  2  J 


a  ? 
p-^r 


1+^/kz1I+  p(ro+ap) 

1+//-5//  (i  +  //)2 

(1+//-5//)2 


+  Po  Co  P 


r.- 


V 


Po 

2 


\ 

- aju 

) 


+ 


(ro+g/j)2 

d+^)2 


+  a 


(3.88) 


As  the  compressibility  of  one  phase  of  soil  is  different  from  another  under  the 
pressure,  the  volume  of  a  particular  soil  phase  cannot  be  explicitly  determined.  For  a 
multi-phase  soil  medium  under  pressure,  either  a  volume  fraction  or  a  weight  fraction 
with  respect  to  the  original  soil  volume  may  be  used  to  determine  the  content  of  each 
phase.  If  the  initial  volume  fractions  of  the  air,  water  and  solid  phases  of  soil  are 
respectively  Aa  ,  Aw  ,  and  As ,  and  A* ,  A* ,  and  A*  under  the  pressure,  and  pa  ,  pw  ,  and 
ps  are  the  initial  densities  of  the  corresponding  phases,  the  following  equations  can  be 


obtained  (Qian  and  Wang  1993): 


4+4  +  4=i 
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.(3.89) 


Po  =AsPs+AWPW+AaPa 


.(3.90) 


4  =  4 


\PQ  J 


.(3.91) 


A,.,  =  A 


w  w 


kJ  Pw-po  )  |  ] 


-k~ 


ps: 


\  r  w  W 


.(3.92) 


A  =  A. 


K(ps-p0)  n 

V  psc; 


,-k; 


.(3.93) 


where  p0  is  the  initial  density  of  the  soil,  P0  are  the  initial  pressures,  ks ,  kw ,  and  ka  are 
the  respective  exponents  of  the  specific  entropy  of  the  solid,  water  and  air  phases,  Cw 
and  Cs  are  the  sound  speeds  in  water  and  solid,  and  P  ,  P  and  P  are  calculated  using 


(3.84).  The  soil  density  under  pressure  p  can  be  expressed  as 


p  =  Po(a;+a;+a;)-1 . (3.94) 

If  the  initial  weight  fractions  of  the  air,  water  and  solid  phases  of  soil  are  respectively  Ra  , 
Rw  ,  and  Rs ,  it  can  be  shown  that 


R  +R  +R  _  PaAa  +Ph.4  +PAs  =  Pp_  =1 


P  0 


Po 


(3.95) 


The  specific  energy  E  and  the  specific  volume  V  of  the  soil  under  pressure  can  be 
expressed  in  terms  of  the  weight  fractions  of  the  three  constituent  phases  as  follows 
(Lovetskii  et  al.  1979): 


E  =  Ra  Ea  +  Rw  Ew  +  Rs  Es 


(3.96) 


V  =  RaVa+RwVw+RsVs 
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(3.97) 


The  values  of  the  EOS  parameters  for  saturated  soil  are  given  in  Table  3-3. 
are  also  valid  for  the  dry  soil,  for  example  A  =  0.68  ,  Aw  =  0  ,  and  Aa  =  0.32 


These  values 
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CHAPTER  FOUR  NUMERICAL  ANALYSIS  AND  COMPARISON  WITH 

TEST  DATA 

4.1  INTRODUCTION 

Since  shock  wave  propagation  in  soils  including  interaction  between  fluid  (air) 
and  solid  (soil  or  structures),  numerical  simulation  of  explosion  in  soils  is  complex. 
Differences  in  characteristics  were  observed  from  detonation  in  two  differing  soil  types: 
dry  sand  and  saturated  sand  (Chapman  et  al.  2006,  Gupta  1999).  How  to  deal  with  soil 
properties  in  the  simulation  of  explosion  is  important  to  obtain  reasonably  good 
simulation  results.  Therefore,  there  are  two  most  important  factors  need  to  be  considered 
for  getting  a  good  simulation.  Two  parameters  are  key  in  dealing  with  soil  properties  in 
explosion  simulation  and  equation  of  state  used. 

Since  the  air  and  water  are  trapped  within  soil  voids  and  deformed  with  the  soil 
skeleton  under  blast  loading,  relative  movement  between  the  skeleton  and  the  water  and 
air  can  be  neglected.  Therefore,  a  stress  tensor  may  be  decomposed  into  a  deviatoric 
stress  component  and  a  hydrostatic  pressure: 

a^a.-PS, . (4.1) 

where  crij  is  the  total  stress,  Pis  hydrostatic  pressure,  positive  in  compression,  anddV  is 
Kronecker  delta.  Deviatoric  stress  can  be  derived  from  soil  material  model  and 


hydrostatic  pressure  can  be  determined  from  an  equation  of  state. 
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Two  methods  are  currently  used  to  consider  the  soil  properties  in  explosion 
simulation,  the  empirical  method  and  the  soil-model  method. 

In  the  empirical  method,  an  equivalent  input  load  is  directly  applied  on  concerned 
structures  while  the  interaction  between  soil  and  explosive/structures  is  neglected.  For 
example,  when  analyzing  a  plate  subjected  to  the  explosion  detonated  from  a  shallow- 
buried  landmine,  an  empirical  relationship  of  a  specific  impulse  (Westine  et  al  1985)  may 
be  directly  applied  on  the  plate;  this  is  known  as  US  Army  TACOM  impulse  model.  The 
main  advantage  of  this  method  is  its  ease  in  application.  Validation  of  this  method  on 
some  simple  geometrical  structures  was  done  with  carefully  calibrated  parameters  in  the 
impulse  model  (Williams  et  al  2002). 

For  the  sake  of  simplification,  the  conventional  way  is  to  apply  an  equivalent 
input  loads  based  on  empirical  functions  which  includes  soil  properties  without  equation 
of  state.  For  example,  *LOAD-BLAST  boundary  condition  was  implemented  into  LS- 
DYNA  finite  element  code  based  on  CONWEP  air-blast  functions  (Randers-Pehrson  and 
Bannister  1997)  to  simulate  surface  detonations.  This  input  load  cannot  consider  the 
effects  on  different  soil  types.  A  more  accurate  empirical  relationship,  called  US  Army 
TACOM  impulse  model,  was  developed  by  Westine  et  al.  (1985)  at  Southwest  Research 
Institute  to  predict  the  impulse  applied  by  a  buried  mine  to  a  plate  at  a  given  offset  from 
the  mine.  The  relationship  is  expressed  as 

in  =  f(r .  d,  D^n  e ,  s,psoil ,  rnmm  t,p,0) . (4.2) 

where  the  soil  density,  psoii,  is  considered.  Other  variables  are  defined  in  FIG.  4-1. 
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FIG.  4-  1  Definition  of  variables  in  US  Army  TACOM  impulse  model 

(Westine  et  al.  1985) 

If  the  model  parameters  are  carefully  calibrated  (Williams  et  al  2002),  this 
empirical  model  can  predict  the  effect  mine  blast  on  simple  geometries  reasonable  well. 

This  method  is  obviously  not  capable  of  capturing  the  complex  transient 
interactions  between  the  soil  and  detonation  products,  which  may  substantially  affect  the 
estimated  blast  loads  and  the  resultant  structure.  Soil  and  debris  could  not  be 
implemented  directly  into  the  soil  finite  element  model. 

In  order  to  compensate  this  limitation,  in  the  soil-model  method  proposed  herein 
the  constitutive  models  are  invoked  to  simulate  the  soil  behavior  in  explosion  (Gupta 
1999,  Wang  2001). 

The  soil  and  foam  material  model  was  applied  by  Wang  (2001)  in  LS-DYNA 
(*MAT5)  to  simulate  a  series  of  explosions  in  air  and  soil.  The  simulation  results  were 
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compared  reasonably  well  with  experimental  results.  The  soil  and  foam  model  may  be 
considered  as  a  special  kind  of  cap  model,  but  the  cap  is  a  plane  cap  in  principal  stress 
space  (Krieg  1972).  Although  this  model  is  highly  efficient,  it  has  the  following 
disadvantages:  lack  of  associative  flow  plasticity,  instability  in  unconfined  states,  and  no 
consideration  of  strain  rate  effect. 

To  date  the  equations  of  state  that  can  be  used  in  numerical  simulation  of 
explosion  in  soils  limited.  Sedgwick  (1974)  applied  for  Tillotson  equation  of  state  in  the 
two-dimensional  HELP  computer  code  to  solve  the  interaction  between  buried  explosive 
charges.  Dynamic  material  properties  experiments  were  performed  to  provide  the 
necessary  soil  equation  of  state  parameters  which  are  required  as  input  to  the  numerical 
model.  The  equation  of  state  for  the  solid  component  and  the  substance  in  the  pores  (gas 
or  liquid)  were  taken  in  the  Mie-gruneisen  form  by  Lovetskill,  Maslennikov  and  Fetisov 
(1979).  The  gaseous  component  was  assumed  to  be  an  ideal  gas  and  the  temperature  of 
all  the  components  was  assumed  to  be  identical.  A  particular  form  of  the  Mie-gruneisen 
equation  of  state  was  applied  by  Grujicic  et  al  (2008)  to  calculate  pressure  dependence  on 
mass  density  and  internal-energy  density.  Qian  (1993)  and  Wang  (2004)  both  applied 
Kandaur  conceptual  equation  of  state  based  on  the  three-phase  soil  structure  in  the  soil 
model  for  blast  loading. 

In  this  chapter,  viscoplastic  cap  soil  model  and  equation  of  state  model  are 
integrated  into  LS-DYNA  finite  element  code  (PC  version)  as  user-defined  material  and 
EOS  model  respectively.  A  series  of  landmine  explosion  tests  in  dry  sand  and  saturated 
sand  conducted  by  Materials  Sciences  Corporation  (2007)  are  simulated  using  the  user- 
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defined  soil  model  and  EOS  model.  The  simulation  results  are  evaluated  through 


comparison  against  experimental  data. 


4.2  PROPERTIES  OF  SOIL  USED  IN  EXPLOSIVE  TESTS 


The  soil  subjected  to  the  plate  impact  tests  by  Chapman  et  al.  (2006)  was  quartz 
sand  provided  by  the  Concrete  Structure  Section  (CSS),  Department  of  Civil  and 
Environmental  Engineering,  Imperial  College,  UK.  The  sand  had  average  particle  size  of 
230  pm  and  dry  soil  density  of  1520  ±  50  kg  m"  .  Since  the  density  of  quartz  is  2650  kg 
m"  ,  the  porosity  of  the  sand  was  about  43%.  If  all  the  voids  were  filled  with  water,  the 
theoretical  maximum  water  content  and  density  would  be  22%  and  1950  kg  m~\ 
respectively. 


A  sandy  soil  was  provided  by  the  Army  Research  Laboratory  (ARL),  Aberdeen, 
MD,  for  the  explosive  tests  conducted  by  Materials  Sciences  Corporation  (2006).  Table 
4-1  provides  a  comparison  of  the  soil  properties.  Since  the  properties  of  the  ARL  soil  are 
very  comparable  to  those  of  the  CSS  quartz  sand,  the  EOS  models  based  on  the  CSS 
quartz  sand  test  data  were  used  in  the  numerical  simulations  of  the  explosive  tests. 


Table  4-  1  Properties  of  soil  specimens 


Soil 

Provided 

by 

Density 
(kg  m'3) 

Volume  ratio 
of  water 

Porosity 

Dry  Sand 

CSS 

1520 

0% 

43% 

ARL 

1871 

0% 

31.23% 

Saturated 

Sand 

CSS 

1950 

22% 

43% 

ARL 

2072 

20.12% 

31.23% 
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4.3  DESCRIPTION  OF  EXPLOSION  TEST 

To  verify  the  validity  of  the  revised  soil  model  under  blast  loading,  the  EOS 
models  along  with  the  viscoplastic  cap  model  are  incorporated  into  the  software  FS- 
DYNA  (LSTC  2003)  as  user-defined  subroutines  for  numerical  simulations. 

Explosive  tests  at  a  3-cm  depth  of  burial  (DOB)  for  dry  (3  tests)  and  saturated  (3 
tests)  sandy  soil  were  conducted  by  the  Materials  Sciences  Corporation  (2006).  Tests  data 
were  provided  by  ARL.  As  shown  in  FIG.  4-2,  a  70-cm  high  cylindrical  tank,  made  of  a 
1.2-cm  thick  steel  pipe  with  a  60-cm  inner  diameter,  was  filled  with  the  test  soil.  A  100- 
gram  C4  explosive  charge  with  6.4-cm  diameter  and  2-cm  thickness  was  placed  at  a  3-cm 
depth  in  the  soil  at  the  center  of  the  tank.  Nine  “pencil”  pressure  transducers  were  placed 
above  the  soil  mass  to  measure  air  pressure  from  the  explosive  gas  bubble  expansion. 
Transducers  #1  through  #5  were  placed  at  the  same  standoff  distance  of  30  cm  and 

O 

pointing  toward  the  center  of  C4  at  0,  ±22.5,  and  ±45  angles,  #6  through  #8  were  placed 

O  O 

at  70  cm  and  at  0  and  ±30  angles,  and  #9  at  1 13  cm  and  at  0  angle.  Transducers  #1,  #6 
and  #9,  respectively  located  at  30  cm,  70  cm  and  110  cm  directly  above  the  soil,  arc 
selected  for  comparisons  between  the  numerical  results  and  measured  air  pressure  due  to 
buried  explosions.  The  scheme  of  the  explosive  tests  set-up  is  shown  in  FIG.  4-3.  FIG.  4- 
4  and  FIG.  4-5  are  explosive  tests  photos  taken  by  high  speed  video  for  saturated  soil  and 


dry  soil,  respectively. 
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FIG.  4-  2  Explosive  test  set-up 


FIG.  4-  3  Schematic  explosive  test  set-up 


40psec  417psec  833psec 


1042psec  1202psec  1455psec 

FIG.  4-  4  Explosive  test  for  saturated  soil  with  DOB=3  cm  by  high  speed  video 


614psec  836psec  1044psec 


FIG.  4-  5  Explosive  test  for  dry  soil  with  DOB=3  cm  by  high  speed  video 


4.4  FINITE  ELEMENT  MODEL 
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Taking  advantage  of  symmetry,  only  a  quarter  of  the  test  setup  was  modeled.  The 
finite  element  model  is  shown  in  FIG.  4-6  containing  a  110-cm  air  volume  above  and  a 
70-cm  soil  volume  below  the  soil  surface,  meshed  with  6,400  8-node  solid  ALE  elements. 
Fine  mesh  was  generated  for  the  explosive  and  for  the  air  and  soil  volumes  surrounding 
the  C4  where  high  strain  gradients  are  anticipated.  The  fine  mesh  of  soil  elements 
extended  3  cm  above  and  below,  and  4.8  cm  outward  in  the  radial  direction  beyond  the 
circumference  of  the  C4  explosive.  The  fine  mesh  of  air  elements  extended  8  cm  above 
the  soil  surface  and  8  cm  in  the  radial  direction.  Coarser  mesh  was  used  in  the  region 
further  away  from  the  explosive  to  reduce  computation  time.  The  materials  used  in  finite 
element  model  and  their  equation  of  states  are  shown  in  FIG.  4-7. 


FIG.  4-  6  Finite  element  mesh 
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Material:  MAT_NULL 

EOS:  EOS  LINEAR-PLOYNOMAL 


Material:  MAT  48  (user  defined) 
EOS:  NO  EOS 


Soil  Material:  MAT  48  (user  defined) 

EOS:  EOS  23  (user  defined) 


C4 


Material:  MAT  HIGH  EXPLOSIVE  BURN 
EOS:  EOS  LINEAR-PLOYNOMAL 


FIG.  4-  7  Material  and  EOS  model 


The  steel  tank  was  treated  as  a  fixed  boundary  of  the  soil.  All  the  exterior 
boundary  of  the  air  was  also  fixed.  The  height  of  the  air  in  the  finite  element  model  was 
set  110cm,  which  was  sufficient  for  investigating  pressure  vs.  time  history  at  the 
positions  of  the  transducers.  The  nodes  on  the  interfaces  between  the  air,  soil  and 
explosive  were  merged,  which  was  the  most  reliable  and  economic  way  to  simulate 
contact. 

To  avoid  large  distortions  in  the  elements  by  the  explosion,  automatic  rezoning 
was  achieved  by  using  the  Arbitrary  Lagrangian-Eulerian  (ALE)  technique  (“LS-DYNA” 
1998).  Set  as  multiple  materials,  explosive,  soil  and  air  were  allowed  within  the  same 
mesh  so  that  the  explosive  product  (i.e.,  the  fire  ball)  would  be  able  to  expand  into  initial 
soil  and  air  meshes  and  the  soil  could  be  ejected  into  the  air  mesh. 

There  are  a  total  of  12  material  parameters  in  the  viscoplastic  cap  model:  rj,  N,  f' 
in  the  viscous  flow  rule;  W,  D.  R,  XQ  in  the  cap  surface;  a,  J3,  /,  0  in  the  failure  surface; 
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and  T  in  the  tension  cutoff  surface.  In  addition,  the  bulk  modulus  K  and  the  shear 
modulus  G  are  needed  for  the  elastic  soil  response.  These  parameters  are  determined 
from  various  static  soil  tests.  Values  of  the  model  parameters  for  a  sandy  soil  are  given 
in  Table  4-2. 


Table  4-  2  Viscoplastic  cap  model  parameters  for  sandy  soil 


Sand 

K  (MPa) 

G  (MPa) 

A  (MPa) 

B  (MPa'1) 

T  (MPa) 

0 

R 

Dry 

106.4 

63.85 

0.0642 

0.34283 

0.00589 

0.18257 

5.00 

Saturated 

1000 

20.00 

0.0625 

0.36430 

0.00320 

0.24900 

5.32 

Sand 

W 

D  (MPa'1) 

Xo(MPa) 

T  (MPa) 

H  (psec'1) 

fo  (MPa) 

N 

Dry 

0.2142 

0.00952 

0.01 

0.0069 

2x1  O’4 

l.OxlO5 

1.0 

Saturated 

0.2250 

0.00884 

0.01 

0.0072 

lxlO’4 

1.2xl05 

1.0 

The  explosion  product  of  C4  is  modeled  with  the  JWL  equation  of  state  (Dobratz 
and  Crawford  1985).  It  can  be  written  in  the  form 


f 

P  =  A  1 

v 


CO 

RyV 


-R,V 


( 

+  B  1 

V 


CO 

R^V 


J 


coE 

V 


(4.3) 


where  A,  B.  Rj,  R2  and  co  are  constants  determined  from  the  experiments,  V  is  the  relative 
volume,  E  is  the  internal  energy.  The  EOS  parameters  for  C4  are  listed  in  Table  4-3. 


Table  4-  3  JWL  equation  of  state  parameters  for  C4 


A  (MPa) 

B  (MPa) 

Ri 

r2 

CO 

E0  (MPa) 

Vo 

609970 

12950 

4.5 

1.4 

0.25 

9000 

1 
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The  air  above  the  soil  is  modeled  with  the  LINEAR-POLYNOMIAL  equation  of 
state  (“LS-DYNA”  1998).  It  can  be  written  in  the  form 


P  —  Cg  +  Cj/y  +  C^/j  +  C3jU~  +  (c4  +  C3ju  +  C6/y  . (4.4) 

where  Co,  Cj,  C2,  C3,  C5  and  Q  are  polynomial  equation  coefficient,  /j  =  —  - 1 ,  and  — 

Pa  Pa 


is  the  ratio  of  current  density  to  reference  density.  E  is  the  internal  energy,  V  is  the 
relative  volume.  The  EOS  parameters  for  air  are  listed  in  Table  4-4. 


Table  4-  4  LINEAR-POLYNOMIAL  equation  of  state  parameters  for  air 


Co 

Q 

c2 

c3 

c4 

c5 

c6 

E0  (MPa) 

Vo 

-1.0e-6 

0.0 

0.0 

0.0 

0.4 

0.4 

0.0 

0.25 

1 

As  illustrated  in  FIG.  4-8,  at  detonation  (time  t  =  0),  energy  prescribed  by  Eq.  (4.3) 
is  released  from  the  center  of  the  C4  elements.  This  pressure  is  transferred  to  the  soil 
elements  surrounding  the  C4,  which  are  within  the  fine  mesh  of  the  model.  The  EOS 
models  developed  are  used  to  account  for  thermodynamic  equilibrium  for  the  air,  water 
and  solid  phases  of  these  soil  elements.  The  shock  front  pressure  decays  exponentially 
with  distance  from  the  point  of  detonation,  and  pressure  level  is  much  lower  when  the 
shock  front  reaches  the  fine  mesh  boundary.  Thus,  EOS  models  are  not  used  for  soil 


elements  in  the  coarse  mesh.  This  process  is  illustrated  in  FIG.  4-9,  simply. 


Air  Elements 


■st  t  t  t  t  \  JM 


FIG.  4-  8  Material  and  EOS  models  used  for  the  FE  mesh 


/ 

[  EOS 


\ 


Interface 


C4 


II 

II 

/ 


Water 


Air 


/ 


Apparent  Properties  -  3  Phases  Combined 


1 


Viscoplastic  Cap  Model  of  Soi 


FIG.  4-  9  Energy  transmission  chart 
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The  flowchart  in  FIG.  4-10  illustrates  the  implementation  of  the  algorithm  using 


an  incremental  time-step  approach.  The  model  is  subjected  to  gravity  load  to  provide  the 
initial  pressure  and  energy  of  the  soil.  The  change  in  volume  over  a  time  step  is 
calculated  for  each  soil  element  after  detonating  the  C4.  The  changes  in  volume  of  the 
three  phases  are  calculated  by  Eq.  (3.91),  (3.92)  and  (3.93).  During  each  time  step,  the 
internal  energy  consisting  of  the  deviatoric  strain  energy  and  the  mechanical  work  done 
by  the  hydrostatic  pressure  is  updated.  The  new  work  done  by  the  pressure  on  the  change 
in  volume  from  each  phase  is  added  to  the  internal  energy  of  the  soil  element  by  Eq. 
(3.96).  The  soil  bulk  modulus  is  updated  using  Eq.  (3.88)  for  subsequent  soil  stress  and 
strain  calculations  in  the  viscoplastic  cap  model.  The  instructions  for  the  implementation 
of  a  user-defined  EOS  are  given  in  the  Appendix  B  of  the  LS-DYNA  user’s  manual 


(LSTC  2003). 
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FIG.  4-  10  Flowchart  showing  the  solution  algorithm  for  use  in  LS-DYNA 
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4.5  SIMULATION  FOR  SATURATED  SOIL 

4.5.1  SIMULATION  CASES  AT  DIFFERENT  ELEMENT 

Case  1  (under  the  C4): 

A  soil  element  (#654),  shown  in  FIG.  4-11,  whose  center  is  located  at  3  cm  to  the 
right  and  2.75  cm  down  from  the  center  of  C4,  is  selected  from  a  saturated  soil  test  to 
illustrate  the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec,  its  soil 
density  is  2055  kg/nr,  bulk  modulus  is  lOOOMPa,  and  the  volume  fractions  of  solid, 
water  and  air  are  respectively  70%,  20%  and  10%.  When  the  shock  arrives  at  time  step  t 
=20+5=25psec,  hydrostatic  pressures  in  the  solid,  water  and  air  phases  are  calculated  to 
be  5.02MPa,  0.0874MPa  and  0.000215MPa,  respectively,  by  Eq.  (3.85),  (3.86)  and  (3.87). 
The  volume  fractions  in  soil  are  updated  using  Eq.  (3.91),  (3.92)  and  (3.93),  to  be  70.15%, 
20.10%  and  9.1%.  Using  Eq.  (3.88)  and  (3.94)  to  update  the  soil  bulk  modulus  and 
density  are  1142.12MPa  and  2063  kg  m"  .  The  soil  volume  increment  can  be  obtained 
from  LS-DYNA,  total  volume  increment  is  -1.7076  E-05  (p=7.601E-05),  solid  volume 
increment  AVS  is  -8.201E-06  (p=4.743E-05),  water  volume  increment  AVW  is  -1.353E-06 
(p=4.10E-05)  and  air  volume  increment  AVa  is  -7.522E-6  (p=3.0852E-03).  The  soil 
pressure  is  2.43MPa.  It  can  be  passed  to  deviatoric  stress  to  calculate  total  stress  by  Eq. 
(4.1).  The  soil  internal  energy  is  0.0000417MPa  by  Eq.  (3.96).  By  now,  all  parameters  of 
viscoplastic  cap  model  and  EOS  are  known.  The  next  time  step  can  be  run.  At  t=  40psec, 
soil  bulk  modulus  arrives  peak  value  2000MPa.  The  volume  fractions  in  soil  are  72.92%, 


20.51%  and  8.75%,  respectively.  The  increments  of  volume  fractions  in  soil  are  2.92%, 
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0.51%  and  -1.25%,  respectively.  The  procedure  of  volume  fractions  change  is  shown  in 


FIG.  4-12. 


FIG.  4-11  Element  #654 
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FIG.  4-  12  Saturated  soil  increments  of  volume  fractions  in  element  #654 


Case  2  (flush  with  the  C4): 


A  soil  element  (#748),  shown  in  FIG.  4-13,  whose  center  is  located  at  3  cm  to  the 
right  from  the  center  of  C4  and  flush  with  the  center  of  C4,  is  selected  from  a  saturated 
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soil  test  to  illustrate  the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec, 
its  soil  density  is  2055  kg/m  ,  bulk  modulus  is  lOOOMPa,  and  the  volume  fractions  of 


solid,  water  and  air  are  respectively  70%,  20%  and  10%.  When  the  shock  arrives  at  time 
step  t  =20+5=25 psec,  hydrostatic  pressures  in  the  solid,  water  and  air  phases  are 
calculated  to  be  3.48MPa,  0.0615MPa  and  0.0002MPa,  respectively,  by  Eq.  (3.85),  (3.86) 
and  (3.87).  The  volume  fractions  in  soil  are  updated  using  Eq.  (3.91),  (3.92)  and  (3.93), 
to  be  70.13%,  20.10%  and  9.08%.  Using  Eq.  (3.88)  and  (3.94)  to  update  the  soil  bulk 
modulus  and  density  are  1 133.23MPa  and  2060  kg  m"  .  The  soil  volume  increment  can  be 
obtained  from  LS-DYNA,  total  volume  increment  is  -1.662  E-05  (p=7.711E-05),  solid 
volume  increment  AVS  is  -8.175E-06  (p=4.743E-05),  water  volume  increment  A Vw  is  - 
1.212E-06  (p=4.10E-05)  and  air  volume  increment  AVa  is  -7.233E-6  (p=3.0852E-03). 
The  soil  pressure  is  2.15MPa.  It  can  be  passed  to  deviatoric  stress  to  calculate  total  stress 
by  Eq.  (4.1).  The  soil  internal  energy  is  0.0000417MPa  by  Eq.  (3.96).  By  now,  all 
parameters  of  viscoplastic  cap  model  and  EOS  are  known.  The  next  time  step  can  be  run. 
At  t=  40psec,  soil  bulk  modulus  arrives  peak  value  2000MPa.  The  volume  fractions  in 
soil  are  72.31%,  20.51%  and  8.84%,  respectively.  The  increments  of  volume  fractions  in 
soil  are  2.31%,  0.51%  and  -1.26%,  respectively.  The  procedure  of  volume  fractions 


change  is  shown  in  FIG.  4-14. 


98 


"" 

X 

X 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/  /  / 

7  7/ 

FIG.  4-13  Element  #748 
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FIG.  4-  14  Saturated  soil  increments  of  volume  fractions  in  element  #748 


Case  3  (above  the  C4): 

A  soil  element  (#852),  shown  in  FIG.  4-15,  whose  center  is  located  at  3  cm  to  the 
right  and  2.75  cm  above  from  the  center  of  C4,  is  selected  from  a  saturated  soil  test  to 
illustrate  the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec,  its  soil 
density  is  2055  kg/m ,  bulk  modulus  is  lOOOMPa,  and  the  volume  fractions  of  solid, 
water  and  air  are  respectively  70%,  20%  and  10%. 
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There  is  a  little  difference  between  the  elements  above  C4  and  the  elements  under 


the  C4.  The  equation  to  calculate  hydrostatic  pressure  can  be  expressed  as: 


P  =  PqC2  h  +  y0Ev . (4.5) 

For  solid,  this  equation  can  be  expressed  as: 

p  =  (2650)(6.319)2//  +  (l.0)£r . (4.6) 

For  water,  this  equation  can  be  expressed  as: 

p  =  (l000)(l.460)2//  +  (0.6)£1’ . (4.7) 

For  air,  this  equation  can  be  expressed  as: 

p  =  (l.2)(0.24l)> . (4.8) 


When  the  shock  arrives  at  time  step  t  =20+5=25psec,  hydrostatic  pressures  in  the 
solid,  water  and  air  phases  are  calculated  to  be  -4.21MPa,  -2.01  IMPa  and  -0.832MPa, 
respectively,  by  Eq.  (4.6),  (4.7)  and  (4.8).  The  volume  fractions  in  soil  are  updated  using 
Eq.  (3.91),  (3.92)  and  (3.93),  to  be  68.99%,  19.44%  and  6.99%.  Using  Eq.  (3.88)  and 
(3.94)  to  update  the  soil  bulk  modulus  and  density  are  923.23MPa  and  2032  kg  m 3.  The 
soil  volume  increment  can  be  obtained  from  LS-DYNA,  total  volume  increment  is  4.233 
E-05  (p=-9.348E-05),  solid  volume  increment  AVS  is  0.875E-05  (p=-3.691E-05),  water 
volume  increment  A Vw  is  1.226E-05  (p=-5.421E-05)  and  air  volume  increment  AVa  is 
2.122E-05  (p=-7.188E-03).  The  soil  pressure  is  -5.786MPa.  It  can  be  passed  to  deviatoric 
stress  to  calculate  total  stress  by  Eq.  (4.1).  The  soil  internal  energy  is  0.0000417MPa  by 
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Eq.  (3.96).  By  now,  all  parameters  of  viscoplastic  cap  model  and  EOS  are  known.  The 
next  time  step  can  be  run.  At  t=  35psec,  since  shock  wave  arrives,  soil  above  C4  is 
blown  by  the  force  of  the  explosion.  The  volume  fractions  in  soil  are  0.0%,  0.0%  and 
0.0%,  respectively.  The  increments  of  volume  fractions  in  soil  are  -70.0%,  -20.0%  and  - 
10.0%,  respectively.  The  procedure  of  volume  fractions  change  is  shown  in  FIG.  4-16. 
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FIG.  4- 15  Element  #852 
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FIG.  4-16  Saturated  soil  increments  of  volume  fractions  in  element  #852 


Case  4  (Air  element  above  the  C4): 
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An  air  element  (#4498),  shown  in  FIG.  4-17,  whose  center  is  located  at  30  cm 
above  from  the  center  of  C4,  is  selected  from  a  saturated  soil  test  to  illustrate  the 
numerical  procedure.  Initially,  this  element  is  defined  by  material  model  of  air.  Since  the 
ALE  (Arbitrary  Lagrangian-Eulerian)  calculation  is  selected  by  this  study,  the  primary 
advantage  of  ALE  is  the  number  and  types  of  materials  present  in  an  element  can  change 
dynamically  when  elements  with  more  than  one  material.  Under  blasting  loading,  a  part 
of  volume  of  the  element  4498  is  occupied  by  soil  debris  following  an  explosion,  shown 
in  FIG.  4-18.  The  soil  volume  fraction  arrive  peak  value  17.1%  at  300psec. 


114498 


FIG.  4-17  Air  element  #4498 
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FIG.  4-  18  Volume  fraction  of  saturated  soil  in  air  element  #4498 


4.5.2  COMPARISON  OF  SIMULATION  WITH  TEST  RESULTS 

From  FIG.  4-19  to  FIG.  4-27  present  the  air  pressure  time-histories  at  three  tests 
respectively,  which  were  recorded  by  the  pencil  gages  (see  Fig.  8)  after  a  C4  charge  was 
detonated  in  saturated  sand  at  a  DOB  =  3  cm  (Materials  Sciences  Corporation  2006).  A 
comparison  between  the  predicted  shock  front  air  pressure  and  the  experimental  data 
obtained  at  distances  of  30  cm,  70  cm,  and  110  cm  directly  above  the  soil  is  shown  in 
FIG.  4-28.  The  difference  between  the  numerical  results  and  the  average  test  data  at  30, 
70  and  110-cm  standoff  distances  are  4.5%,  12.5%  and  7.2%,  respectively. 

Density  and  bulk  modulus  are  the  most  sensitive  parameters  in  simulation  model. 
A  comparison  among  simulation  results  with  the  density  decreased  to  90%  of  initial 
density  and  with  the  bulk  modulus  decreased  to  90%  of  initial  bulk  modulus  and  density 
and  bulk  modulus  keep  the  initial  value  is  shown  in  FIG.  4-29. 
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FIG.  4-19  Saturated  sand  air  pressure  time-histories,  30  cm  standoff  distance  #1 
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FIG.  4-  20  Saturated  sand  air  pressure  time-histories,  70  cm  standoff  distance  #1 
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FIG.  4-21  Saturated  sand  air  pressure  time-histories,  110  cm  standoff  distance  #1 
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FIG.  4-  22  Saturated  sand  air  pressure  time-histories,  30  cm  standoff  distance  #2 
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FIG.  4-  23  Saturated  sand  air  pressure  time-histories,  70  cm  standoff  distance  #2 
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FIG.  4-  24  Saturated  sand  air  pressure  time-histories,  110  cm  standoff  distance  #2 
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FIG.  4-  25  Saturated  sand  air  pressure  time-histories,  30  cm  standoff  distance  #3 


FIG.  4-  26  Saturated  sand  air  pressure  time-histories,  70  cm  standoff  distance  #3 
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FIG.  4-  27  Saturated  sand  air  pressure  time-histories,  110  cm  standoff  distance  #3 
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28  Comparison  between  numerically  predicted  and  experimental  values  for 
saturated  sand  (Shock  front  pressure  in  air  VS.  Transducer  distance) 
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FIG.  4-  29 


Comparison  of  simulation  results  due  to  parameters  change  for  saturated  soil 
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The  saturated  soil  volume  fractions  of  three  phases  in  a  part  of  finite  element 
mash  before  the  shack  wave  arriving  is  shown  in  FIG.  4-30.  The  saturated  soil  volume 
fractions  of  three  phases  in  a  part  of  finite  element  mash  at  the  180psec  is  shown  in  And 
FIG.  4-31. 


FIG.  4-  30  Saturated  soil  volume  fractions  of  three  phases  before  the  shack  wave 


arrives 


FIG.  4-31  Saturated  soil  volume  fractions  of  three  phases  at  1  BOpsec 

4.6  SIMULATION  FOR  DRY  SOIL 

4.6.1  SIMULATION  CASES  AT  DIFFERENT  ELEMENT 

Case  1  (under  the  C4): 

A  soil  element  (#654),  shown  in  FIG.  4-11,  whose  center  is  located  at  3  cm  to  the 
right  and  2.75  cm  down  from  the  center  of  C4,  is  selected  from  a  dry  soil  test  to  illustrate 
the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec,  its  soil  density  is 
1802  kg/m  ,  bulk  modulus  is  106.4MPa,  and  the  volume  fractions  of  solid,  water  and  air 
are  respectively  68%,  0.0%  and  32%.  When  the  shock  arrives  at  time  step  t 
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=20+5=25  psec,  hydrostatic  pressures  in  the  solid,  water  and  air  phases  are  calculated  to 
be  2.13MPa,  O.OMPa  and  0.000215MPa,  respectively,  by  Eq.  (3.85),  (3.86)  and  (3.87). 
The  volume  fractions  in  soil  are  updated  using  Eq.  (3.91),  (3.92)  and  (3.93),  to  be  70.04%, 
0.0%  and  29.97%.  Using  Eq.  (3.88)  and  (3.94)  to  update  the  soil  bulk  modulus  and 
density  are  117.16MPa  and  1811  kg  m"  .  The  soil  volume  increment  can  be  obtained  from 
LS-DYNA,  total  volume  increment  is  -2.387E-05  (p=4.481E-06),  solid  volume 
increment  AVS  is  -7.879E-06  (p=4.743E-05),  water  volume  increment  AVW  is  0.0  and  air 
volume  increment  AVa  is  -1.599E-05  (p=6.163E-03).  The  soil  pressure  is  1.62MPa.  It  can 
be  passed  to  deviatoric  stress  to  calculate  total  stress  by  Eq.  (4.1).  The  soil  internal  energy 
is  0.0000175MPa  by  Eq.  (3.96).  By  now,  all  parameters  of  viscoplastic  cap  model  and 
EOS  are  known.  The  next  time  step  can  be  mn.  At  t=  40psec,  soil  bulk  modulus  arrives 
peak  value  513MPa.  The  volume  fractions  in  soil  are  82.68%,  0.0%  and  17.4%, 
respectively.  The  increments  of  volume  fractions  in  soil  are  14.68%,  0.0%  and  -14.6%, 
respectively.  The  procedure  of  volume  fractions  change  is  shown  in  FIG.  4-32. 
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FIG.  4-  32  Dry  soil  volume  fraction  in  element  #654 
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Case  2  (flush  with  the  C4): 

A  soil  element  (#748),  shown  in  FIG.  4-13,  whose  center  is  located  at  3  cm  to  the 
right  from  the  center  of  C4  and  flush  with  the  center  of  C4,  is  selected  from  a  dry  soil  test 
to  illustrate  the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec,  its  soil 
density  is  1802  kg/m ,  bulk  modulus  is  106.4MPa,  and  the  volume  fractions  of  solid, 
water  and  air  are  respectively  68%,  0.0%  and  32%.  When  the  shock  arrives  at  time  step  t 
=20+5=25 psec,  hydrostatic  pressures  in  the  solid,  water  and  air  phases  are  calculated  to 
be  3.72MPa,  O.OMPa  and  0.000201MPa,  respectively,  by  Eq.  (3.85),  (3.86)  and  (3.87). 
The  volume  fractions  in  soil  are  updated  using  Eq.  (3.91),  (3.92)  and  (3.93),  to  be  70.02%, 
0.0%  and  30.0%.  Using  Eq.  (3.88)  and  (3.94)  to  update  the  soil  bulk  modulus  and  density 
are  112.34MPa  and  1808  kg  m"  .  The  soil  volume  increment  can  be  obtained  from  LS- 
DYNA,  total  volume  increment  is  -1.933  E-05  (p=6.412E-05),  solid  volume  increment 
AVS  is  -8.119E-06  (p=4.919E-05),  water  volume  increment  A Vw  is  0.0  and  air  volume 
increment  AVa  is  -1.1131E-05  (p=7.075E-03).  The  soil  pressure  is  1.33MPa.  It  can  be 
passed  to  deviatoric  stress  to  calculate  total  stress  by  Eq.  (4.1).  The  soil  internal  energy  is 
0.0000175MPa  by  Eq.  (3.96).  By  now,  all  parameters  of  viscoplastic  cap  model  and  EOS 
are  known.  The  next  time  step  can  be  run.  At  t=  40psec,  soil  bulk  modulus  arrives  peak 
value  513MPa.  The  volume  fractions  in  soil  are  81.33%,  0.0%  and  18.4%,  respectively. 
The  increments  of  volume  fractions  in  soil  are  13.33%,  0.0%  and  -13.6%,  respectively. 


The  procedure  of  volume  fractions  change  is  shown  in  FIG.  4-33. 
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FIG.  4-  33  Dry  soil  volume  fraction  in  element  #748 

Case  3  (above  the  C4): 

A  soil  element  (#852),  shown  in  FIG.  4-15,  whose  center  is  located  at  3  cm  to  the 
right  and  2.75  cm  above  from  the  center  of  C4,  is  selected  from  a  dry  soil  test  to  illustrate 
the  numerical  procedure.  Before  the  shock  wave  arrives  at  t  =20psec,  its  soil  density  is 
1802  kg/m  ,  bulk  modulus  is  106.4MPa,  and  the  volume  fractions  of  solid,  water  and  air 
are  respectively  68%,  0.0%  and  32%.  When  the  shock  arrives  at  time  step  t 
=20+5=25 psec,  hydrostatic  pressures  in  the  solid,  water  and  air  phases  are  calculated  to 
be  -3.62MPa,  O.OMPa  and  0.0313MPa,  respectively,  by  Eq.  (4.6),  (4.7)  and  (4.8).  The 
volume  fractions  in  soil  are  updated  using  Eq.  (3.91),  (3.92)  and  (3.93),  to  be  65.03%,  0.0% 
and  30.10%.  Using  Eq.  (3.88)  and  (3.94)  to  update  the  soil  bulk  modulus  and  density  are 
103.02MPa  and  1800  kg  m"  .  The  soil  volume  increment  can  be  obtained  from  LS- 
DYNA,  total  volume  increment  is  5.119  E-05  (p=-7.762E-05),  solid  volume  increment 
AVS  is  1.496E-05  (p=-4.743E-05),  water  volume  increment  AVW  is  0.0  and  air  volume 
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increment  AVa  is  3.623E-05  (p=-4.0852E-03).  The  soil  pressure  is  -0.94MPa.  It  can  be 
passed  to  deviatoric  stress  to  calculate  total  stress  by  Eq.  (4.1).  The  soil  internal  energy  is 
0.0000175MPa  by  Eq.  (3.96).  By  now,  all  parameters  of  viscoplastic  cap  model  and  EOS 
are  known.  The  next  time  step  can  be  run.  At  t=  35psec,  since  shock  wave  arrives,  soil 
above  C4  is  blown  by  the  force  of  the  explosion.  The  volume  fractions  in  soil  are  0.0%, 
0.0%  and  0.0%,  respectively.  The  increments  of  volume  fractions  in  soil  are  -68.0%,  0.0% 
and  -32.0%,  respectively.  The  procedure  of  volume  fractions  change  is  shown  in  FIG.  4- 
34. 


Time  (usee) 

FIG.  4-  34  Dry  soil  volume  fraction  in  element  #852 


Case  4  (Air  element  above  the  C4): 

An  air  element  (#4498),  shown  in  FIG.  4-17,  whose  center  is  located  at  30  cm 
above  from  the  center  of  C4,  is  selected  from  a  dry  soil  test  to  illustrate  the  numerical 
procedure.  Under  blasting  loading,  a  part  of  volume  of  the  element  4498  is  occupied  by 
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soil  debris  following  an  explosion,  shown  in  FIG.  4-35.  The  soil  volume  fraction  arrive 
peak  value  2.17%  at  250psec. 


Time  (psec  E+03) 

FIG.  4-  35  Volume  fraction  of  dry  soil  in  air  element  #4498 


4.6.2  COMPARISON  OF  SIMULATION  WITH  TEST  RESULTS 

From  FIG.  4-36  to  FIG.  4-44  present  the  air  pressure  time-histories,  which  were 
recorded  by  the  pencil  gages  (see  Fig.  8)  after  a  C4  charge  was  detonated  in  dry  sand  at  a 
DOB  =  3  cm  (Materials  Sciences  Corporation  2006). 

The  predicted  shock  front  arrival  times  in  the  air  directly  above  the  explosion  are 
compared  against  those  read  from  the  data  traces  recorded  by  transducers  #1,  #6  and  #9 
in  FIG.  4-45.  The  difference  between  the  predicted  shock  front  arrival  time  and  the 


average  test  data  at  0,  22.5  and  45  offset  angles  are  1.8%,  4.4%,  and  9.7%,  respectively. 
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A  comparison  between  the  predicted  shock  front  air  pressure  and  the  experimental 
data  obtained  at  distances  of  30  cm,  70  cm,  and  110  cm  directly  above  the  soil  is  shown 
in  FIG.  4-46.  The  difference  between  the  numerical  results  and  the  average  test  data  at 
30,  70  and  110-cm  standoff  distances  are  2.2%,  20%  and  64%,  respectively.  A 
comparison  among  simulation  results  when  the  density  is  decreased  to  90%  of  initial 
density,  when  the  bulk  modulus  is  decreased  to  90%  of  initial  bulk  modulus  and  density 
and  bulk  modulus  keep  the  initial  value  is  shown  in  FIG.  4-47. 

The  dry  soil  volume  fractions  of  three  phases  in  a  part  of  finite  element  mash 
before  the  shack  wave  arriving  is  shown  in  FIG.  4-48.  The  dry  soil  volume  fractions  of 
three  phases  in  a  part  of  finite  element  mash  at  the  120psec  are  shown  in  FIG.  4-49. 

The  soil  ejecta  heights  between  high  speed  video  and  numerical  simulation  at  time 
=  420,  830  and  1040  ps  since  detonation  for  tests  in  dry  sand  and  in  saturated  sand  are 
compared  in  FIG.  4-50,  4-51  and  4-52,  respectively.  The  maximum  difference  between 
the  predicted  and  measured  ejecta  heights  is  24%  for  explosive  tests  in  dry  sand  and  9.6% 
in  saturated  sand. 


FIG.  4-  36  Dry  sand  air  pressure  time-histories,  30  cm  standoff  distance  #1 
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FIG.  4-  37  Dry  sand  air  pressure  time-histories,  70  cm  standoff  distance  #1 
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FIG.  4-  38  Dry  sand  air  pressure  time-histories,  1 10  cm  standoff  distance  #1 


FIG.  4-  39  Dry  sand  air  pressure  time-histories,  30  cm  standoff  distance  #2 
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FIG.  4-  40  Dry  sand  air  pressure  time-histories,  70  cm  standoff  distance  #2 
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FIG.  4-  41  Dry  sand  air  pressure  time-histories,  1 10  cm  standoff  distance  #2 


FIG.  4-  42  Dry  sand  air  pressure  time-histories,  30  cm  standoff  distance  #3 
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FIG.  4-  43  Dry  sand  air  pressure  time-histories,  70  cm  standoff  distance  #3 
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FIG.  4-  44  Dry  sand  air  pressure  time-histories,  1 10  cm  standoff  distance  #3 
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FIG.  4-  45  Comparison  between  numerically  predicted  and  experimental  values  for  dry 


sand  (Blast  wave  arrival  time  VS.  Transducer  offset  angle) 
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FIG.  4-  46  Comparison  between  numerically  predicted  and  experimental  values  for  dry 
sand  (Shock  front  pressure  in  air  VS.  Transducer  distance) 
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FIG.  4-  47  Comparison  of  simulation  results  due  to  parameters  change  for  dry  soil 
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FIG.  4-  48  Dry  soil  volume  fractions  of  three  phases  before  the  shack  wave  arrives 
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FIG.  4-  49  Dry  soil  volume  fractions  of  three  phases  at  120psec 
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Dry  Sand  D0B=3  cm 
Time=420  nsec 


Saturated  Sand  D0B=3  cm 
Time  =417  nsec 


FIG.  4-  50  Comparison  of  soil  ejecta  heights:  High  speed  video  vs.  Simulation 


at  time  =  420  psec 
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Dry  Sand  D0B=3  cm 
Time=S36  psec 


Saturated  Sand  D0B=3  cm 
Time  =S33  psec 


FIG.  4-51  Comparison  of  soil  ejecta  heights:  High  speed  video  vs.  Simulation 

at  time  =  830  psec 
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Saturated  Sand  D0B=3  cm 
Time  =1042  psec 


FIG.  4-  52  Comparison  of  soil  ejecta  heights:  High  speed  video  vs.  Simulation 

at  time  =  1040  psec 


4.7  CONCLUSIONS 
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By  means  of  comparison  against  experimental  data,  predictions  of  the  viscoplastic 
cap  model  demonstrate  better  agreement  than  those  of  the  inviscid  cap  model,  more 
accurate  7%  than  inviscid  cap  model,  since  the  viscoplastic  model  can  capture  the  high 
strain-rate  (with  durations  in  milliseconds)  effects  on  explosion  simulation.  The  high 
strain  effects  are  manifested  by  an  apparent  increase  of  shock  wave  propagation  speed, 
peak  overpressure  and  impulse.  Although  the  effects  on  certain  variables  are  not  apparent, 
such  as  the  air  blast  shock  wave  propagation  and  the  explosion  characteristics,  the  high 
strain  rate  effects  are  generally  significant  (Jackson  et  al,  1980)  and  cannot  be  neglected 
in  explosion  simulation. 

A  fine  mesh  (about  0.14cm  per  element)  needs  to  be  used  in  order  to  improve  the 
simulation  accuracy.  Besides,  the  high  strain-rate  effects  need  to  be  studied  through 
explosion  tests  in  clayey  soils  in  order  to  draw  a  general  conclusion. 
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CHAPTER  FIVE  CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  CONCLUSIONS 

This  thesis  investigates  and  proposes  soil  models,  implements  the  models  in  LS- 
DYNA  finite  element  code,  and  evaluates  the  performance  of  viscoplastic  cap  material 
model  with  equation  of  states  through  comparison  against  available  explosion  test  data. 
The  soil  behavior  under  blast  loading,  the  phenomena  due  to  explosion,  are  particularly 
studied. 

Two  viscoplastic  cap  models  are  based  on  Perzyna’s  theory  and  Duvant-Lions’ 
theory  studied.  By  comparing  with  the  solutions  to  a  hypothetical  loading  test,  the  two 
viscoplastic  models  produce  virtually  identical  responses  when  the  viscous  parameters 
are  judiciously  selected  for  each  model.  However,  differences  between  the  Perzyna’s  and 
the  Duvant-Lions’  model  were  observed  when  simulating  the  experiment  tests  conducted 
under  rapid  loading.  The  prediction  of  the  Perzyna’s  model  appears  to  agree  better  (about 
4%)  with  experimental  data  than  that  of  the  Duvant-Lions’  model,  and  the  Perzyna’s 
model  is  more  flexible  for  data  fitting,  more  accurate  about  6.6%  than  Duvant-Lions’ 
model.  Therefore,  the  Perzyna’s  viscoplastic  cap  model  is  implemented  into  LS-DYNA 
to  represent  the  soil  model  with  consideration  of  strain-rate  effect. 

To  improve  the  accuracy  of  simulation  results,  three  phases  equation  of  states  are 
developed  based  on  Mie-Gruneison  equation  of  state.  For  soil  mass  surrounding  the 
source  of  energy  release,  equation  of  state  models  for  the  three  phases  of  soil  are 


developed  as  each  of  the  three  phases  responds  differently  to  shock  loading.  Finally, 
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these  three  phase  equation  of  states  have  been  integrated  with  the  viscoplastic  cap  model 
and  incorporated  into  the  LS-DYNA  software  as  user-supplied  subroutines  for  numerical 
simulations  of  explosive  tests  in  dry  soil  as  well  as  in  wet  soil. 

By  means  of  comparison  against  experimental  data,  the  predicted  time  of  arrival 
and  the  overpressure  in  air  directly  above  buried  explosions  agree  well  with  the 
experimental  data.  There  was  noticeable  improvement  using  the  revised  cap  model  with 
EOS  for  the  prediction  of  wet  soil  behavior  under  blast  loading  than  dry  soil.  It  is 
concluded  that  the  revised  cap  model  with  EOS  is  adequate  for  blast  loading  behavior 
simulations  for  soil  with  different  degrees  of  water  contents. 

5.2  RECOMMENDATIONS 

Refinements  of  viscoplastic  cap  models  would  include  a  more  realistic  treatment 
for  tension  cutoff  for  the  Perzyna  type  and  a  more  elaborate  formulation  for  the  Duvant- 
Lions  type.  The  former  is  very  important  to  the  simulation  of  underground  explosion.  The 
latter  is  to  improve  the  flexibility  of  the  Duvant-Lions’  model. 

A  series  of  explosion  tests  needs  to  be  conducted  in  clayey  soil  with  different 
degrees  of  water  contents.  As  evidenced  by  the  previous  experimental  studies,  the  clayey 
soils  are  more  sensitive  to  the  loading  rate  than  the  sandy  soil.  If  these  tests  are  being 
conducted,  the  comprehensive  static  tests  for  the  same  soil  should  also  be  conducted  to 
calibrate  soil  model  parameters  and  EOS  parameters.  More  simulations  can  be  ran  with 
equation  of  state  for  soils  with  various  degrees  of  saturation.  This  step  is  essential  for 
ensuring  the  accuracy  of  the  numerical  simulations. 
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SUBROUTINE  OF  USER  DEFINED  MATERIAL  MODEL 


subroutine  umat48(cm,eps,sig,epsp,hsv,dtl,capa,etype,tt, temper, failel,crv) 
c  Perzyna's  Viscoplastic  Cap  Model  for  Soil : 
c  cm(l)  =  young's  modulus 

c  cm(2)  =  possion's  ratio 

c  cm(3)  =  buckling  modulus 

c  cm(4)  =  shear  modulus 

c  cm(5)  =  alfa  in  Faliure  Surface 

c  cm(6)  =  beta  ... 

c  cm(7)  =  gama  ... 

c  cm(8)  =  theta  ... 

c  cm(9)  =  r  cap  surface  axis  ratio 

c  cm(10)=  d  hardening  law  exponent 

c  cm(l  1)=  w  hardeng  law  coefficient  (limit  plastic  strain) 
c  cm(12)=  xO  initial  hardening  pressure 

c  cm(13)=  tcut  tension  cut  off  (tcut<=0) 

c  cm(14)=  conv  convegent  factor  (default  value  =  0.001) 

c  cm(15)=  itmat  maximum  iteration  number  (default  value  =  1000) 
c  hsv(l)=total  z-component  strain 

c  hsv(2)=hardening  parameter,  kn 

c  hsv(3)=volumetric  plastic  strain,  evpn 

c  hsv(4)=first  stress  invariant,  J 1 

c  hsv(5)=second  deviatoric  stress  invariant,  SJ2 

c  hsv(6)=response  mode  number,  mode 

include  'iounits.inc' 
character* (*)  etype 

dimension  cm(*),eps(*),sig(*),hsv(*),crv(101,2,*) 

dimension  cmat(6,6),dmat(6,6),hr(6,6),hh(6,6),dfds(6),ddfdds(6,6) 

dimension  ddfdsl(6),dfaids(6),ab(6),sigl(6),se(6) 

real*4  kn,knl,ln,lnl,knl0,k0 

logical  faille 

c  Input  the  user  defined  material  parameters 
bulk=cm(l) 
gshr=cm(2) 
alfa=cm(3) 
beta=cm(4) 
gama=cm(5) 
theta=cm(6) 


r=cm(7) 

d=cm(8) 

w=cm(9) 

x0=cm(10) 

tcut=cm(ll) 

conv=cm(12) 

itmax=cm(13) 

yita=cm(14) 

faiO=cm(15) 

expon=cm(16) 

c  Calculate  the  initial  hardening  parameter  kO  or  input  kn 
if  (hsv(2).eq.O)  then 

if  (xO.ge.  10000.0)  then 
kn=xO 
else 

call  capi(xO,r,alfa, beta, gama, theta, kO,ieer) 
if  (ieer.eq.10)  then 
kO=xO 
endif 
kn=kO 
endif 

else 

kn=hsv(2) 

endif 

c  Form  the  elastic  material  matrix  [cmat]  and  its  reverse  matrix  [dmat] 
cmatii=bulk+4.0/3.0*gshr 
cmatij=bulk-2.0/3.0*gshr 
cmatii=gshr 
do  140  1=1,6 
do  140  j=l,6 

140cmat(i,j)=0.0 
do  160  1=1,6 

if  (i.le.3)  then 
cmat(i,i)=cmatii 
do  150  j=l,3 

if  (i.ne.j)  cmat(i,j)=cmatij 
150  continue 
else 

cmat(i,i)=cmatjj 
endif 
160  continue 

call  mrevs(6, cmat, dmat) 


c  Calculate  the  elastic  trial  strss  { sigl }  =  { sigO}  +  [cmat]: {eps} 
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do  180 1=1,6 
aa=0.0 
do  170  j=l,6 
aa=aa+cmat(i,j)*eps(j) 

170  continue 

sig  1  (i)=-(sig(i)+aa) 

180  continue 

c  Given  other  initial  values 
if  (kn.ge.  10000.0)  then 
xn=x0 
else 

xn=kn+r*(alfa-gama*exp(-beta*kn)+theta*kn) 

endif 

evpn0=w*(l-exp(-d*(xn-x0))) 

dlamd=0.0 

dk=0.0 

knl=kn 

c  Deal  with  tension  cutoff 

sj  1  e=sig  1  ( 1  )+sig  1  (2)+sig  1  (3) 
if  (sjle.gt.(-tcut))  goto  450 
sjlnl=-tcut 
ppe=sjl  e/3.0 
ppt=sjlnl/3.0 
dse=0.0 
do  190  i=l,6 

if  (i.le.3)  then 
se(i)=sigl(i)-ppe 
fmu=1.0 
else 

se(i)=sigl(i) 

fmu=2.0 

endif 

dse=dse+fmu*se(i)*se(i) 

190  continue 

sj2e=sqrt(0.5*dse) 

sj2nl=sj2e 

sj2t=alfa-gama*exp(-beta*(-tcut))+theta*(-tcut) 
if  (sj2e.gt.sj2t)  sj2nl=sj2t 
ratio=0.0 

if  (sj2e.ne.0.0)  ratio=sj2nl/sj2e 
do  200  i=  1,6 
if  (i.le.3)  then 
sig  1  (i)=se(i)*ratio+ppt 
else 


sigl(i)=se(i)*ratio 
endif 
200  continue 
goto  800 

c  Check  other  status  of  the  elastic  trial  stress 
450  continue 

call  differ(sigl,knl,fai,dfds,ddfdds,ddfdsl,dfaids, dfaidl, hsk 
$  , mode, alfa, beta, gama, theta, r,d,w,x0,tcut,yita,fai0,expon) 

if  (mode.eq.0)  goto  800 
residi=fai-dlamd/yita/dt  1 

c  ***  local  iteration  to  fulfill:  residi  =  fai  -  dlamd/yita/dt  =>  convergence 
itt=0 

500  itt=itt+l 

c  2.1:  [hh]  =  I  [cmat]  +  dlamd*  [ddfdds]  1-1 
do  510  i=l,6 
do  510  j=l,6 

510  hr(i,j)=dmat(i,j)+dlamd*ddfdds(i,j) 
call  mrevs(6,hr,hh) 

c  2.2:  divd  =  { dfaids }:  [hh] :{ dlamd*  {ddfdsl}+{dfds} }  +  1/yita/tt  -  dfaidl 
divd=0.0 

do  520  i=l,6 
ab(i)=0.0 
do  520  j=l,6 

ab(i)=ab(i)+hh(i,j)*(dlamd*ddfdsl(j)+dfds(j)) 

520  continue 
do  530  i=l,6 

divd=divd+ab(i)*dfaids(i) 

530  continue 

divd=divd+ 1 .0/yita/dt  1  -dfaidl 
c  2.3:  dlamd  =  dlamd  +  residi/divd  ; 
dlamd=dlamd+residi/divd 

c  2.4:  { sigl }  =  { sig}  +  [cmat]*{ { eps} -dlamd* {dfds} } 
devpn=0.0 
do  550  i=l,6 
ac=0.0 
do  540  j= 1,6 

ac=ac+cmat(i,j)*(-eps(j)-dlamd*dfds(j)) 

540  continue 

sigl(i)=-sig(i)+ac 

if  (i.le.3)  devpn=devpn+dlamd*dfds(i) 

550  continue 


if  (kn.ge.  10000.0)  then 
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knl=kn 

else 

evpn 1 =evpnO+devpn 
if  (evpnl.ge.w)  evpnl=0.9*w 
xn  1 =- 1 .0/d*log(  1 .O-evpn  l/w)+xO 
knlO=knl 
itk=0 

570  continue 
itk=itk+l 

ff=kn  1 0+r*(alfa-gama*exp(-beta*kn  1 0)+theta*kn  1 0)-xn  1 
if  (abs(ff).lt.abs(conv*knlO))  goto  580 
dfekn=gama*beta*exp(-beta*kn  1 0)+theta 
kn  1  =kn  1 0-ff/(  1 .0+r*dfekn) 
knl0=knl 

if  (itk.ge.itmax)  then 

c  write(6,*)'not  convergence  for  knl-kn-ff,knl,kn,ff 

goto  580 

endif 

goto  570 
580  continue 
endif 

c  2.5:  residi(new)  =  fai  -  dlamd/yita/dt 

call  differ(sig  1  ,kn  1  ,fai,dfds,ddfdds,ddfdsl,dfaids,dfaidl,hsk 
$  , mode, alfa, beta, gama, theta, r,d,w,x0,tcut,yita,fai0,expon) 

residi=fai-dlamd/yita/dt  1 

c  2.6:  check  if  the  convergence  condition  is  satisfied 
if  (abs(residi).lt.abs(conv))  goto  800 
if  (itt.ge.itmax)  then 
write(6,*)  'NOT  CONVERGE', itt,residi 
else 

goto  500 
endif 

c  ***  local  iteration  end 
800  continue 
c  Output  the  variables 
do  820  i=l,6 
820  sig(i)=-sigl(i) 


hsv(  1  )=hsv(  1  )+eps(3) 
hsv(2)=knl 

c  write(6,*)'output',mode,knl,kn 

xn  1  =kn  1  +r*(alfa-gama*exp(-beta*kn  1  )+theta*kn  1 ) 
hsv(3)=hsv(3)+w*(1.0-exp(-d*(xnl-x0))) 
pn  1 =(sig(  1  )+sig(2)+sig(3  ))/3 .0 
hsv(4)=pn  1*3.0 


dsigll=sig(l)-pnl 
dsig22=sig(2)-pn  1 
dsig33=sig(3)-pnl 

dsa=dsigl  l*dsigl  l+dsig22*dsig22+dsig33*dsig33 
dse=sig(4)**2+sig(5)**2+sig(6)**2 

hsv(5)=sqrt(0.5*dsa+dse) 

hsv(6)=float(mode) 

return 

end 

c  determine  the  initial  hardening  parameter  kO  according  to  xO 
c 

subroutine  capi_dup(xO,r,alfa, beta, gama, theta, kO,ieer) 
real*4  kn,kO 
ieer=0 

cretia=le-5*(alfa-gama) 

itc=0 

k0=0.0 

40  fO=alfa-gama*exp(-beta*kO)+theta*kO 
dfekO=gama*beta*exp(-beta*kO)+theta 
f=k0+r*f0-x0 
if  (abs(f).lt.cretia)  goto  60 
kn=k0-f/(  1 .0+r*dfek0) 
k0=kn 
itc=itc+l 

if  (itc.gt.60)  goto  50 
goto  40 
50  ieer=10 
60  return 
end 

c 

c  calcuate  the  reverse  matrix 
c 

subroutine  mrevs(ns,cmat,dmat) 

dimension  cmat(ns,ns)  ,dmat(ns,ns) 
do  100  i=l,ns 
do  100j=l,ns 
100  dmat(i,j)=cmat(i,j) 
do  200  n=l,ns 
diag=dmat(n,n) 
do  130  j=l,ns 

130  dmat(n,j)=-dmat(n,j)/diag 
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do  150  i=l,ns 

if  (n.eq.i)  goto  150 
do  140j=l,ns 
if  (n.eq.j)  goto  140 

dmat(i,j)=dmat(i,j)+dmat(i,n)*dmat(n,j) 

140  continue 
150  dmat(i,n)=dmat(i,n)/diag 
dmat(n,n)= 1 .0/diag 
200  continue 
return 
end 

c  calcuate  the  flow  vector  of  failure  surface 

subroutine  differ(ssn  1  ,kn  1  ,fai,dfds,ddfdds,ddfdsl,dfaids,dfaidl, 
$  hsk, mode, alfa, beta, gama, theta, r,d,w,x0,tcut,yita,fai0,expon) 
real*4  knl,lnl 

dimension  ssnl(6),dfds(6),ddfdds(6,6),ddfdsl(6),dfaids(6) 
dimension  dj  lds(6),dj2ds(6),se(6),amat(6,6) 
c  Get  the  basic  flow  vector  :  dj  Ids  =  { dj  1/ds } ;  dj2ds  =  { dj2/ds } 
sj  1 =ssn  1  ( 1  )+ssn  1  (2)+ssn  1  (3) 
pp0=sj  1/3.0 
toth=2. 0/3.0 
aa=0.0 
do  120  i=l,6 
if  (i.le.3)  then 
djlds(i)=1.0 
else 

djlds(i)=0.0 

endif 

se(i)=ssnl(i)-dj  lds(i)*pp0 
if  (i.le.3)  then 
dj2ds(i)=se(i) 
else 

dj2ds(i)=2.0*se(i) 

endif 

aa=aa+se(i)  *di  2ds(i) 
do  100  j=l,6 
amat(i,j)=0.0 
if  (i.eq.j)  then 
if  (i.le.3)  amat(i,j)=toth 
if  (i.gt.3)  amat(i,j)=2.0 
else 

if  (i.le.3.and.j.le.3)  amat(i,j)=-0.5*toth 
endif 
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100  continue 
120  continue 
sj2=sqrt(0.5*aa) 

c  Check  the  stress  status:  Mode  =  0  ->  elastic  ;  -1  ->  tension; 
c  1  ->  failure  ;  2  ->  cap 

mode=0 
fval=0.0 
dfj  1=0.0 
dfj2=0.0 
ddfj  1=0.0 
ddfj2=0.0 
ddfj  12=0.0 
ddfj  lk=0.0 
ddfj2k=0.0 
dkdl=0.0 
dfdk=0.0 
In  1  =max(kn  1,0.0) 
c  if  (sjl.le.-tcut)  then 
c  write(6,*) 'differ- l',sj  l,-tcut,lnl 

c  fval=si2-(-tcut) 

c  dfj  1=1.0 

c  if  (fval.gt.0.0)  mode=-l 

c  else  if  (sjl.gt.-tcut.and.sjl.le.ini)  then 

if  (sjl.le.lnl)  then 

fval=sj2-(alfa-gama*exp(-beta*sj  l)+theta*sj  1) 
dfj  l=-gama*beta*exp(-beta*sj  l)-theta 
dfj2=0.5/sj2 

ddfj  l=gama*beta*beta*exp(-beta*sj  1) 
ddfj  2=-0 . 25/sj  2/sj  2/sj  2 
if  (fval.gt.0.0)  mode=l 
else 

xn  1  =kn  1  +r*(alfa-gama*exp(-beta*kn  1  )+theta*kn  1 ) 

al=(sjl-lnl)/r 

alr=al/r 

a2=(xnl-lnl)/r 

aa=sqrt(a  1  *a  1 +sj  2*  sj  2) 

a3=1.0/aa/aa/aa 

fval=aa-a2 

dfjl=al/aa/r 

dfj2=0.5/aa 

ddfj  1 =-  a  1  r  *  a  1  r  *  a3 + 1  .0/aa/r/r 
ddfj2=-0.25*a3 
ddfjl2=-0.5*alr*a3 
dldk=0.0 

if  (knl.gt.O)  dldk=1.0 
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ddfj  lk=-ddfj  1  *dldk 
ddfj2k=-ddfj  12*dldk 
dfedl=gama*beta*exp(-beta*kn  1  )+theta 
dkdl=3.0*dfjl/(w*d*exp(-d*(xnl-x0)))/(1.0+r*dfedl) 
dfdk=-dfj  1  *dldk-dfedl*dldk 
if  (fval.gt.O.O)  mode=2 
endif 

if  (fval.le.O.O)  goto  800 
c 

c  MODE  !=  0  — >  viscoplasticity 

c  dfai=dfai/df  ;  dfaidl=dfai/dlamd  ;  dfds=df/ds  ;  ddfdds=ddf/dds 
c 

fai=(fval/faiO)  *  *expon 
dfai=expon*(fval/faiO)**(expon- 1 .0)/fai0 
hsk=dkdl 

dfaidl=dfai*dfdk*dkdl 
do  140  i=l,6 

dfds(i)=dfjl*dj  Ids(i)+dfj2*dj2ds(i) 
ddfdsl(i)=(ddfj  lk*dj  Ids(i)+ddfj2k*dj2ds(i))*dkdl 
dfaids(i)=dfai*dfds(i) 

140  continue 

do  160  i=l,6 
do  160  j=l,6 

ddfdds(i,j)=ddfj  1  *dj  lds(i)*dj  lds(j)+ddfj  12*(dj  lds(i)*dj2ds(j)+ 
$  dj  Ids(j)*dj2ds(i))+ddfj2*dj2ds(i)*dj2ds(j)+dfj2*amat(i,j) 

160  continue 
800  return 
End 
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SUBROUTINE  OF  USER  DEFINED  EOS  MODEL 


subroutine  ueos23s(iflag,cb,pnew,hist,rho0,eosp,specen, 
&  df,dvol,v0, pc, dt,tt,crv, first) 

c***  variables 


c  iflag - =0  calculate  bulk  modulus 

c  =1  update  pressure  and  energy 

c  cb - bulk  modulus 

c  pnew - new  pressure 

c  hist - history  vai’iables 

c  rhoO - reference  density 

c  eosp - EOS  constants 

c  specen  —  energy/reference  volume 

c  df - volume  ratio,  v/vO  =  rhoO/rho 

c  dvol - change  in  volume  over  time  step 

c  vO - reference  volume 

c  pc - pressure  cut-off 

c  dt - time  step  size 

c  tt - current  time 

c  crv - curve  array 

c  first - logical  .true,  for  tt, crv, first  time  step 

c  (for  initialization  of  the  history  variables) 

c 


include  'nlqparm' 
logical  first 

dimension  hist(*),eosp(*),crv(101,2,*) 
real*4  As,Aw,Aa,dvols,dvolw,dvola 

c  solid, water, air— precent 

AsO  =eosp(l) 

AwO  =eosp(2) 

AaO  =eosp(3) 

c  solid, water, air— density 

rs  =eosp(4) 
rw  =eosp(5) 
ra  =eosp(6) 

c  solid, water, air— ks,kw,ka 

sk  =eosp(7) 
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wk  =eosp(8) 
ak  =eosp(9) 

c  input  parameters— mixed  soil 
c  =eosp(10) 
si  =eosp(ll) 
s2  =eosp(12) 
s3  =eosp(13) 
gO  =eosp(14) 
sa  =eosp(15) 
si  1 =s  1 - 1 . 
s22=2.*s2 
s33=3.*s3 
s32=2.*s3 
sad2=.5*sa 
g0d2=l.-.5*g0 

roc2=rho0*c**2 

c  input  parameters— solid 
cs  =eosp(16) 
ssl  =eosp(17) 
ss2  =eosp(18) 
ss3  =eosp(19) 
gsO  =eosp(20) 
ssa  =eosp(21) 
ssll=ssl-l. 
ss22=2.*ss2 
ss33=3.*ss3 
ss32=2.*ss3 
sads2=.5*ssa 
g0ds2=l.-.5*gs0 
rocs2=rs*cs**2 

c  input  pai’ameters— water 
cw  =eosp(22) 
swl  =eosp(23) 
sw2  =eosp(24) 
sw3  =eosp(25) 
gwO  =eosp(26) 
swa  =eosp(27) 
swll=swl-l. 
sw22=2.*sw2 
sw33=3.*sw3 
sw32=2.*sw3 
sadw2=.5*swa 


gOdw2=l.-.5*gwO 

rocw2=rw*cw**2 

input  parameters- -air 

ca  =eosp(28) 

sal  =eosp(29) 

sa2  =eosp(30) 

sa3  =eosp(31) 

gaO  =eosp(32) 

saa  =eosp(33) 

sall=sal-l. 

sa22=2.*sa2 

sa33=3.*sa3 

sa32=2.*sa3 

sada2=.5*saa 

gOda2=l.-.5*gaO 

roca2=ra*ca**2 

pO  =eosp(34) 

if  (hist(l).eq.O)  then 
hist(l)=AsO 
hist(2)=AwO 
hist(3)=AaO 
hist(4)=0.0 
hist(5)=0.0 
hist(6)=0.0 
As=hist(l) 
Aw=hist(2) 
Aa=hist(3) 
dvols=hist(4) 
dvolw=hist(5) 
dvola=hist(6) 
else 

As=hist(l) 

Aw=hist(2) 

Aa=hist(3) 

dvols=hist(4) 

dvolw=hist(5) 

dvola=hist(6) 

endif 

RRs=As*rs/rhoO 

RRw=Aw*i-w/rhoO 

RRa=Aa*ra/rhoO 


specens=specen*RRs 

specenw=specen*RRw 

specena=specen*RRa 


Vsold=(df*v0-2.0*dvol)*As 
V  wold=(df  *v0-2  .O*dvol)  *  Aw 
Vaold=(df*v0-2.0*dvol)*Aa 

***  calculate  the  bulk  modulus  for  the  EOS  contribution  to  the  sound  speed 
if  (iflag.eq.O)  then 
from  solid 
xmu=1.0/df-l. 
dfmu=df*xmu 
facp=.5*(l.+sign(l.,xmu)) 
facn=l.-facp 

xnum=l.+xmu*(+g0d2-sad2*xmu) 

xdem=l.-xmu*(sll+dfmu*(s2+s3*dfmu)) 

tmp=facp/(xdem*xdem) 

a=roc2*xmu*(facn+tmp*xnum) 

b=gO+sa*xmu 

pnum=roc2*(facn+facp*(xnum+xmu*(g0d2-sa*xmu))) 
pden=2.*xdem*(-sl  1  +dfmu*(-s22+dfmu*(s2-s33+s32*dfmu))) 
cb=pnum*(facn+tmp)-tmp*a*pden+sa*specen+ 

&  b*df**2*max(pc,(a+b*specen)) 

if  (cb.lt. 0.02)  then 
cb=cb 
else 

cb=0.02 

endif 

***  update  the  pressure  and  internal  energy 
else 

from  solid 
dfs=df*(As/AsO) 
xmus=1.0/dfs-l. 
dfmus=dfs*xmus 
facps=.5*(l.+sign(l.,xmus)) 
facns=l.-facps 

xnums=l.+xmus*(+g0ds2-sads2*xmus) 

xdems=l.-xmus*(ssll+dfmus*(ss2+ss3*dfmus)) 

tmps=facps/(xdems*xdems) 

a=rocs2*xmus*(facns+tmps*xnums) 


b=gsO+ssa*xmus 

dvovOs=.5*(dvols)/(AsO*vO) 

denoms=l.+b*dvovOs 

pnews=(a+specens*b)/max(l.e-6,denoms) 

pnews=max(pnews,pc) 

specens=specens-pnews*dvovOs 


from  water 
dfw=df *  ( A  w/A  wO) 
xmuw=  1 .0/dfw- 1 . 
dfmuw=dfw*xmuw 
facpw=.5*(l.+sign(l.,xmuw)) 
facnw=l.-facpw 

xnumw=l.+xmuw*(+g0dw2-sadw2*xmuw) 

xdemw=l.-xmuw*(swll+dfmuw*(sw2+sw3*dfmuw)) 

tmpw=facpw/(xdemw*xdemw) 

a=rocw2*xmuw*(facnw+tmpw*xnumw) 

b=gwO+swa*xmuw 

dvovOw= .  5  *  (dvolw)/(  A  wO  *  vO) 

denomw=l.+b*dvovOw 

pneww=(a+specenw*b)/max(l.e-6,denomw) 

pneww=max(pneww,pc) 

specenw=specenw-pneww*dvovOw 

from  air 

dfa=df*  (Aa/AaO) 

xmua=1.0/dfa-l. 

dfmua=dfa*xmua 

facpa=.5*(l.+sign(l.,xmua)) 

facna=l.-facpa 

xnuma=  1  .+xmua*(+g0da2-  sada2  *xmua) 

xdema=l.-xmua*(sall+dfmua*(sa2+sa3*dfmua)) 

tmpa=facpa/(xdema*xdema) 

a=roca2*xmua*(facna+tmpa*xnuma) 

b=gaO+saa*xmua 

dvovOa= .  5  *  (dvola)/(  AaO  *  vO) 

denoma=l  .+b*dvovOa 

pnewa=(a+specena*b)/max(  1  .e-6,  denoma) 

pnewa=max(pnewa,pc) 


specena=specena-pnewa*dvovOa 


if  (pnews/=0.0.AND.pneww/=0.0.AND.pnewa/=0.0.AND.dvol/=0.0)  then 
pnew=(pnews*dvols+pneww*dvolw+pnewa*dvola)/dvol 
else 

pnew=pnews+pneww+pnewa 

endif 

specen=specens+specenw+specena 

As=As*(sk*(pnews-pO)/(rs*cs**2)+l)**(-(sk)**(-l)) 

Aw=Aw*(wk*(pneww-pO)/(rw*cw**2)+l)**(-(wk)**(-l)) 

Aa=Aa*(pnewa/pO)**(-(ak)**(-l)) 

hist(l)=As 

hist(2)=Aw 

hist(3)=Aa 

dvols=df*vO*As-Vsold 

dvolw=df*vO*Aw-Vwold 

dvola=df*vO*Aa-Vaold 

hist(4)=dvols 

hist(5)=dvolw 

hist(6)=dvola 


endif 

return 

end 
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LS-DYNA  INPUT  DECK 


$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

$  LS—DYNA (971)  DECK  WRITTEN  BY  :  eta/FEMB-PC  version  28.0 
$  ENGINEER  : 

$  PROJECT  : 

$  UNITS  :  MM,  TON,  SEC,  N 
$  TIME  :  03:01:01  PM 
$  DATE  :  FRIDAY,  Aug  20,  2010 

$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

* KEYWORD 

$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

*TITLE 


dob32.dyn  (explosive) 

$  $ 

$  CONTROL  CARD  $ 

$  $ 

*CONTROL_ACCURACY 
$  OSU  INN  PIDOSU 

1  3 

*  CONTROL_ALE 


$ 

DCT 

NADV 

METH 

AFAC 

BFAC 

CFAC 

DFAC 

EFAC 

2 

1 

2 

-1.0 

0.0 

0.0 

0.0 

0.0 

$ 

START 

END 

AAFAC 

VFACT 

PRIT 

EBC 

PREF 

NSIDEBC 

0.0 

o 

o 

o 

o 

o 

o 

o 

o 

0 

o 

o 

•k 

CONTROL_ENERGY 

$ 

HGEN 

RWEN 

SLNTEN 

RYLEN 

2 

2 

1 

1 

*  CONTROL_HOURGLAS S 
$  I HQ  QH 

4  0.1 

*CONTROL_CONTACT 


$  SLSFAC  RWPNAL 

ISLCHK 

SHLTHK 

PENOPT 

THKCHG 

ORIEN 

ENMASS 

1  0.0 

1 

0 

1 

0 

1 

0 

$  USRSTR  USRFRC 

NSBCS 

INTERM 

XPENE 

SSTHK 

ECDT 

TIEDPRJ 

0  0 

10 

0 

4.0 

0 

0 

0 

*C0NTR0L_0UTPUT 

$  NPOPT  NEECHO 

NREFUP 

IACCOP 

OPIFS 

IPNINT 

IKEDIT 

IFLUSH 

0  0 
$  IPRTF 

0 

*CONTROL_TIMESTEP 

0 

0 

o 

o 

I SDO 

TSLIMT 

DT2MS 

LCTM 

ERODE 

MS1ST 

$  DTINIT  TSSFAC 

1.0e-04  0.2 

0 

o 

o 

o 

o 

0 

1 

0 

$  DT2MSF 

0.0 

*C0NTR0L_TERMINATI0N 

$  ENDTIM  ENDCYC 

DTMIN 

ENDENG 

ENDMAS 

1000.0  0 

1.0 

0 . 0 

o 

o 

1 

C\l 

1 

1 

1 

1 

+ 

1 

1 

1 

1 

\ — 1 

1 

1 

1 

1 

+ 

1 

1 

1 

</> 

1 

1 

CO 

1 

1 

1 

1 

+ 

1 

1 

1 

— + - 4- 

1 

1 

1 

+ 

1 

1 

1 

1 

(_n 

1 

- + - 6- 

1 

r- 

l 

l 

l 

l 

+ 

l 

l 

l 

1 

1 

1 

+ 

1 

1 

1 

1 

CO 
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$  $ 

$  DATABASE  CONTROL  FOR  BINARY  $ 

$  $ 

$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 


*  DAT ABAS  E_B I NARY_D 3  P  LOT 

$  DT/CYCL  LCDT  BEAM  NPLTC 

5.0 

*  DAT ABAS  E_B I NARY_D 3THDT 

$  DT/CYCL  LCDT 

5.0 


$ 

$  DATABASE  EXTENT  CARDS 

$ 

*DATABASE_EXTENT_B INARY 

$A 

$  NEIPH  NEIPS  MAXINT  STRFLG  SIGFLG  EPSFLG  RLTFLG  ENGFLG 

15  1 

$  CMPFLG  IEVERP  BEAMIP  DCOMP  SHGE  STSSZ  N3THDT  IALEMAT 

0 

$  NINTSLD 

$  1 

$ 

$  PART  CARDS 

$ 

$ - + - 1 - + 2 - + - 3 - + 4 - + - 5 - + - 6 - + 7 - + 

*PART 

$HEADING 

SOIL 

$  PID  SECID  MID  EOSID  HGID  GRAV  ADPOPT  TMID 

1  1  12  6 

*PART 

$HEADING 

C4 

$  PID  SECID  MID  EOSID  HGID  GRAV  ADPOPT  TMID 

2  2  2  2 

*PART 

$HEADING 


AIR 

$  PID  SECID  MID  EOSID  HGID  GRAV  ADPOPT  TMID 

3  3  3  1 

$  $ 

$  SECTION  CARDS  $ 

$  $ 


*  SECTION _ SOLID _ ALE 


$ 

SECID 

ELFORM 

AET 

1 

11 

$ 

AFAC 

BFAC 

CFAC 

DFAC 

START 

END 

AAFAC 

o 

o 

o 

o 

o 

o 

o 

o 

0.0 

o 

o 

o 

o 

•k 

SECTION_SOLID 

_ALE 

$ 

SECID 

ELFORM 

AET 

2 

11 

$ 

AFAC 

BFAC 

CFAC 

DFAC 

START 

END 

AAFAC 

o 

o 

o 

o 

o 

o 

0.0 

0.0 

o 

o 

0.0 

k 

SECTION_SOLID 

_ALE 

$ 

SECID 

ELFORM 

AET 
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$ 

3 

AFAC 

11 

BFAC 

CFAC 

DFAC 

START 

END 

AAFAC 

o 

o 

o 

o 

o 

o 

o 

o 

0.0 

o 

o 

o 

o 

$ - + - ! - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

$  $ 


$ 

MATERIAL 

CARDS 

$ 

$ 

$  1  1 

-+ - 2- 

_F0AM 

1 

1 

Cl 

1 

1 

1 

1 

+ 

1 

1 

1 

1  4 

l  5 

I  6 — 

$ 

_ | _ Q 

*MAT_SOIL_AND 
$AM— 1 

i  O 

$  MID 

RO 

G 

BULK 

A0 

A1 

A2 

PC 

1 

1 . 8 

0.0006385 

0.303. 

4000E-137 

.  0330E-07 

0.30-6. 

900E-08 

$  VCR 

REF 

0.0 

o 

o 

$  EPS1 

EPS2 

EPS3 

EPS4 

EPS5 

EPS6 

EPS7 

EPS8 

o 

o 

-0 . 104 

-0.161 

-0.192 

-0.224 

-0.246 

-0.271 

-0.283 

$  EPS9 

EPS10 

-0.29 

-0.40 

$  PI 

P2 

P3 

P4 

P5 

P6 

P7 

P8 

o 

o 

0.00020 

0.00040 

0.00060 

0 . 0012 

0 . 0020 

0 . 0040 

0 . 0060 

$  P9 

P10 

0.0080 

0.041 

*MAT_USER_DEFINED_MATERIAL_MODELS 


$  MID 

RO 

MT 

LMC 

NHV 

IORTHO 

IBULK 

IG 

11 

1 . 800000 

46 

16 

10 

0 

3 

4 

$  IVECT 

IFAIL 

I THERM 

IHYPER 

IEOS 

0 

0 

0 

0 

0 

$  E 

MU 

BULK 

G 

ALFA 

BETA 

GAMA 

THETA 

0.0 

O 

o 

0.0010646 

. 3850e-04  6 

. 4200e-05 

3428.30005 

. 8900e-06 

0 . 182500 

$  R 

D 

W 

XO 

TCUT 

CONV 

ITMAX 

5.000000 

952.00000 

0.214200 

0.06 

. 9000e-08 

0 . 001000 

60 . 000000 

0 . 00e+0 

*MAT_USER_DEFINED_MATERIAL_MODELS 

$  MID 

RO 

MT 

LMC 

NHV 

IORTHO 

IBULK 

IG 

12 

2.05500 

48 

16 

9 

0 

1 

2 

$  IVECT 

IFAIL 

I THERM 

IHYPER 

IEOS 

0 

0 

0 

0 

1 

$  BULK 

G 

ALFA 

BETA 

GAMA 

THETA 

R 

D 

0.0100002 

. 0000e-04  6 . 

.  2500e-07 

3643.00003 

. 2000e-08 

0.249000 

5.320000 

0.00884 

$  W 

XO 

TCUT 

CONV 

UMAX 

YITA 

FAI0 

EXPON 

0.225000 

0 . 001e-31 . 

.  2000e-08 

0.001000 

60 . 000000 

0 . 200e-l 

1.2 

1.0 

*MAT _ HIGH _ EXPLOSIVE _ BURN 

$  AM-2 

$  MID 

RO 

D 

PCJ 

BETA 

K 

G 

SIGY 

2 

1 . 601 

0.8193 

0.28 

0.0 

0.0 

0.0 

o 

o 

* INITIAL_DETONATION 

$  PID 

X 

Y 

z 

LT 

2 

0.0 

0.0 

o 

1 

*MAT_NULL 

$AM-3 

$  MID 

RO 

PC 

MU 

TEROD 

CEROD 

YM 

PR 

3 

0.00129 

o 

o 

0.0 

0.0 

0.0 

0.0 

0.0 

*E0S_JWL 

$  AEQUATI0N_2 
$  EOSID 

A 

B 

R1 

R2 

OMEGA 

E0 

VO 

2 

6.0997 

0.1295 

4.5 

1 . 4 

0.25 

0.090 

1.0 

*EOS  LINEAR  POLYNOMIAL 
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$  AEQUATION _ 1 


$  EOSID 

CO 

Cl 

C2 

C3 

C4 

C5 

C6 

1-0. 

0000010 

0.0 

0.0 

0.0 

0 . 40 

0 . 40 

0.0 

$  E0 

vo 

0.0000025 

1.0 

*  E  0  S_GRUNE I S  EN 

$  AEQUATION_4 

$  EOSID 

c 

SI 

S2 

S3 

GAMAO 

A 

E0 

4 

0.032 

4  .  92 

0.0 

0.0 

1 . 11 

0.0 

0.0 

$  VO 

1.0 

*EOS_USER_DEFINED 

$  AEQUATION _ 5 

$  EOSID 

EOST 

LMC 

NHV 

IVECT 

EO 

VO 

BULK 

5 

21 

6 

6 

0 

0.0 

1.0 

0 . 002064 

$  C 

SI 

S2 

S3 

GAMAO 

A 

0.032 

4 . 92 

0.0 

0.0 

0 . 11 

0.0 

*EOS_USER_DEFINED 

$  AEQUATION _ 6 

$  EOSID 

EOST 

LMC 

NHV 

IVECT 

EO 

vo 

BULK 

6 

23 

34 

6 

0 

0.0 

1.0 

0.0 

$  AsO 

AwO 

AaO 

Rs 

Rw 

Ra 

ks 

kw 

0.7 

0.2 

0 . 1 

2 . 65 

1.0 

0 . 0012 

3.0 

7.0 

$  ka 

C 

SI 

S2 

S3 

GAMAO 

A 

C-s 

1 . 4 

0.032 

4  .  92 

0.0 

0.0 

0 . 11 

0.0 

0 . 6319 

$  Sl-s 

S2-s 

S3-s 

GAMAO-s 

A-s 

C-w 

Sl-w 

S2-w 

1.41 

0.0 

0.0 

1.0 

0.0 

0.146 

2.0 

0.0 

$  S3-w 

GAMAO-w 

A-w 

C-a 

Sl-a 

S2-a 

S3-a 

GAMAO-a 

0.0 

0.6 

0.0 

0.02406 

1 .0602 

0.0 

0.0 

0.4 

$  A-a 

PO 

0.0 

t _ | _ _ 

1.0e-07 

$ 

_ 1 _ 2  — 

_ I _ 2  — 

_ | _ ^ _ 

_ l _ 2  — 

_ | _ g  _ 

_ I _ '7_ 

_ | _ g 

$ 

$ 

SEGMENT  SET 

CARDS 

$ 

$  $ 
$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 


*  SET_SEGMENT 
$  A  S  E  GMENT_S  E  T  1 
$  SID 

1 

DAI 

0.0 

DA2 

0.0 

DA3 

0.0 

DA4 

0.0 

$  N1 

N2 

N3 

N4 

A1 

A2 

A3 

A4 

3278 

3292 

7092 

7086 

o 

o 

o 

o 

o 

o 

o 

o 

6870 

6876 

6900 

6894 

o 

o 

o 

o 

o 

o 

o 

o 

7816 

7817 

7799 

7798 

o 

o 

o 

o 

o 

o 

o 

o 

7817 

7818 

7800 

7799 

o 

o 

o 

o 

o 

o 

o 

o 

+ 

1 

1 

1 

1 

\ — 1 

1 

1 

1 

1 

+ 

1 

1 

1 

<j>  <j>  <j> 

- 2- 

1 

cn 

i 

1 

1 

1 

+ 

1 

i 

i 

NODE  SET 

CARDS 

1 

1 

1 

+ 

1 

1 

1 

1 

Ol 

1 

1 

1 

1 

+ 

1 

1 

1 

1 

<r>  <r>  co 

$  $ 
$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 


*SET_ 

t  A 

_N0DE _ LIST 

9 

$ 

SID 

DAI 

DA2 

DA3 

DA4 

1 

o 

o 

o 

o 

o 

o 

o 

o 

$ 

NIDI 

NID2 

NID3 

NID4 

NID5 

NID6 

NID7 

NID8 

121 

122 

123 

124 

129 

130 

135 

136 

141 

142 

145 

148 

153 

154 

157 

160 

160 


1665 

1667 

1675 

1677 

1679 

1681 

1683 

1685 

7651 

7667 

7668 

7672 

7689 

7691 

7693 

7722 

7724 

7726 

7728 

7730 

7746 

7747 

7751 

7768 

7770 

7772 

7801 

7803 

7805 

7807 

7809 

*SET_ 

NODE_LIST 

$ASPC 

CARD  AT 

NODE  SET 

2 

$ 

SID 

DAI 

DA2 

DA3 

DA4 

2 

o 

o 

o 

o 

o 

o 

o 

o 

$ 

NIDI 

NID2 

NID3 

NID4 

NID5 

NID6 

NID7 

NID8 

13 

14 

15 

16 

23 

24 

31 

32 

39 

40 

44 

48 

55 

56 

60 

64 

71 

72 

76 

80 

87 

88 

92 

96 

7758 

7760 

7777 

7779 

7781 

7783 

7785 

*SET_ 

NODE_LIST 

$ASPC 

CARD  AT 

NODE  SET 

3 

$ 

SID 

DAI 

DA2 

DA3 

DA4 

3 

o 

o 

o 

o 

o 

o 

o 

o 

$ 

NIDI 

NID2 

NID3 

NID4 

NID5 

NID6 

NID7 

NID8 

213 

214 

223 

232 

241 

250 

259 

3854 

3863 

3872 

472 

473 

482 

491 

500 

509 

518 

1261 

1277 

1293 

1309 

1325 

1341 

1740 

1756 

1772 

1788 

1804 

1820 

3561 

3577 

4049 

4058 

4067 

4076 

4085 

4094 

4383 

4392 

4401 

4410 

4419 

4428 

4437 

4446 

4455 

4464 

5143 

5159 

5175 

5191 

5207 

5223 

5624 

5640 

5656 

5672 

7755 

5688 

5704 

7360 

7376 

7392 

7597 

7676 

*SET_ 

NODE_LIST 

$ASPC 

CARD  AT 

NODE  SET 

4 

$ 

SID 

DAI 

DA2 

DA3 

DA4 

4 

o 

o 

o 

o 

o 

o 

o 

o 

$ 

NIDI 

NID2 

NID3 

NID4 

NID5 

NID6 

NID7 

NID8 

2195 

2920 

2921 

2945 

2969 

2993 

3017 

3041 

3065 

3089 

3113 

3137 

3161 

3185 

3209 

3233 

t  i 

3580 

3579 

+ - 2 - 

3585 

3583 

3581 

3578 

9 - h 

$ 

_ _ 

- h 

_ g _ 

__| _ ^  _ 

_ I _ 3 _ |_  _ 

_ g  _ 

_ I _ 7  _ 

_ | _ g 

$ 

$ 

BOUNDARY 

NON  REFLECTING  CARDS 

$ 

$ 

$ 

$ - + 

- 1 - 

+ - 2 - 

—  + 

- 3 - 

1 

■5J* 

1 

1 

1 

1 

+ 

1 

i 

1 

i 

+ 

i 

i 

i 

1 

(_n 

1 

1 

1 

1 

+ 

- 6- 

- + - 7- 

- + - 8 

*BOUNDARY_NON_REFLECTING 
$  ANON-REFLECTING  CARD  1 
$  SSID  AD  AS 

1  0.0  0.0 


$ - + - ! - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

$  $ 

$  BOUNDARY  SPC  CARDS  $ 

$  $ 

$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

*  BOUND  ARY_S  P  C_S  ET_ID 
$  ID 

1 


$ — + 

$ 

$ 

$ 

$ — + 
*ALE_ 
$AALE 
$ 

*ale_: 

$AALE 

$ 

*ale_: 

$AALE 

$ 

$ - + 

$ 

$ 

$ 

$ - + 

*NODE 

$ 


NSID 

1 

2 

2 

3 

3 

4 
4 


CID  DOFX  DOFY  DOFZ 

0  0  10 


DOFRX 

1 


DOFRY 

0 
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DOFRZ 

1 


MULT I - 
_1 

PSID 

1 

MULT I - 
_2 

PSID 

2 

MULT I - 
_3 
PSID 
3 


0  111 

ALE  CARDS 

-MATERIAL_GROUP 

IDTYPE 

1 

-MATERIAL_GROUP 

IDTYPE 

1 

-MATERIAL_GROUP 

IDTYPE 

1 

NODE  INFORMATION 


1 

+ - 7 


1 - + - £ 


—  7 - + - £ 


1 

1 

1 

+ 

1  X 

1 

1 

1 

CM 

1 

1 

1 

-3 - + - 4- 

Y 

1 

1 

1 

+ 

1 

1 

1 

1 

CJi 

1 

1 

1 

1 

+ 

(SI  1 

1 

1 

-6 - +- 

TC 

1 

1 

1 

<1 

&  1 

n  i 

2.262741 
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$  ELEMENTS  INFORMATION 

$ 

$ 

$  SOLID  ELEMENTS 

$ 

*ELEMENT_SOLID 

$  EID  PID  NIDI  NID2  NID3  NID4  NID5 

1  1  1  2  3  4  5 

2  1  4  3  9  10  8 


-6 - + - 7 - + - £ 

-6 - + - 7 - + - £ 


NID6  NID7  NID8 

6  7  8 

7  11  12 


6398  3  7737  7738  7720  7719  7816  7817  7799  7798 

6399  3  7738  7739  7721  7720  7817  7818  7800  7799 

$ - + - 1 - + - 2 - + - 3 - + - 4 - + - 5 - + - 6 - + - 7 - + - 8 

*END 


